Phase field (PF) method is incorporated into the coupled hydromechanical (HM) model in this paper to simulate the preferential gas flow in saturated bentonite. In order to account for the effect of the heterogeneous distribution of HM properties on the development of preferential pathways, a spatially auto-correlated random field is generated and applied to the coupled HM-PF model. A stabilized technique for equal low-order finite element is adopted in this paper to avoid the pore pressure oscillation and to increase the computational efficiency. The simulated results showed that the developed coupled HM-PF model is able to explicitly describe the preferential gas flow in saturated bentonite. Meanwhile, some of the commonly observed behaviors in experiments, such as the nearly saturated state after gas injection, the heterogeneous distribution of water pressure and the build-up of both water pressure and total stress during fracturing process, are successfully captured by the developed model. In addition, four factors, including the Gaussian random field, the coefficient of variation, the boundary condition and the injection rate, have been examined in terms of their effects on the fracturing process. The numerical results showed that these four factors have significant influences on both the fracture pattern and the evolution of gas pressure.