As a safe long-term management of nuclear wastes, deep geological repositories (DGRs) have been proposed or currently being constructed in several countries. The host rocks in DGRs are saturated with water after the geological disposal facilities (GDFs) are closed and sealed. Significant gas can be generated due to several processes, e.g., the metal corrosion, water radiolysis or microbial reaction of organic materials, etc. The generated gas is anticipated to span throughout the long-term disposal of waste, which may jeopardize the stability of host rocks. Correspondingly, the performance of GDF will be affected since the host rocks provide a final impediment to the radionuclide transport. As gas migration in saturated host rocks is a highly coupled hydro-mechanical (HM) process, either gas-induced micro-fracturing or macro-fracturing may contribute to the development of preferential gas pathways, which needs to be concerned to ensure the feasibility and safety of geological disposal. Current numerical studies on the gas migration behavior devoted to explaining the experimental phenomena in the gas injection tests conducted on the rock materials, in which some behaviors still cannot be well represented, i.e., gas induced fracturing, volulme dilation, anisotropic radial deformation. Therefore, to better represent the actual physical process of preferential gas flow, two modelling frameworks, i.e., macroscopic HM framework and two-scale HM framework, are proposed in the PhD study. For the macroscopic HM framework, a double porosity model is firstly developed based on the dual continuum method, in which the volumetric strains of the porous continuum (PC) and fractured continuum (FC) are work-conjugated to the respective effective stress level. The treatment in two types of porosity allows us to capture that the opening/closure of the fractures is caused by the interaction between the dilation of the PC and the dilation of the FPM, which is beneficial to describe the gas induced fracturing in an implicit way. Then, an enriched embedded fracture model (EFM) is proposed to address the mechanical behavior of fractures. A hyperbolic relation of fracture deformability is incorporated into the rock matrix, as a result the fractured rock shows a nonlinear elastic behavior, which can capture the stiffness degradation due to fracture opening. The equivalent continuum method is provided to derive the effective compliance tensor, which includes the transverse isotropic matrix and two fracture sets. Using the enriched EFM with a three-dimensional (3D) geometry is able to capture the anisotropic radial deformation during gas migration. Although the macroscopic HM framework is able to capture the major HM behaviors related to preferential gas flow, the development of gas dilatant pathways is still represented in an implicit way. Therefore, a two-scale HM framework is developed to explicitly simulate the development of preferential gas pathways. Initiating from the periodically distributed microstructures with microcracks, the asymptotic homogenization method is used to derive the macroscopic governing equations coupled with the normalized damage variable. The time-dependent damage evolution law is obtained from the microscopic mechanical energy analysis for evolving microcracks. Both time effect and size effect are incorporated in the damage model that will affect the overall HM behavior of rocks. The developed two-scale HM framework with single gas flow can qualitatively capture important behaviors, such as the discrete pathways, localized gas flow, unstabilized fracture branching. More specifically, the simulated results demonstrates that the inter-connection of fractures from gas inlet to outlet is a prerequisite for gas breakthrough, accompanied by large amounts of gas flowing out of the sample and a rapid drop in gas injection pressure. Incorporating water flow in the two-scale framework allows the model to quantitatively reproduce the experimental phenomena observed in the laboratory air injection tests, such as gas pressure evolution and mechanical deformation. More importantly, the model exlpaines that the significant differences in controlling gas breakthrough and mechanical deformation are resulting from the arbitrary nature of microstructural heterogeneities. To account for the gas-water interaction in the two-scale HM framework, a fully coupled two-phase flow and elaso-damage model is developed to simulate the laboratory and in-situ gas injection experiments. The model can quantitatively capture the experimental behaviors, e.g., gas pressure evolution and non-desaturation phenomenon. Furthermore, model results show that the highly localized fracture pathways are the major places where gas and water interacts each other, and as a result the rock is still kept fully saturated. As a whole, the obtained numerical results are synthesized and analyzed, the pros and cons of the developed models are discussed. To better improve the model performance, some recommendations are proposed for the future studies.