A thermodynamically consistent phase field model that accounts for initial stress state is proposed in this paper to simulate the gas migration process in saturated bentonite. The energy contribution due to the fracturing process is included in Coussy's thermodynamic framework for unsaturated porous media. The possible effect of the interfaces between different phases on the driving force functional for phase field and the effective stress has been identified from the proposed thermodynamic framework. In addition, the initial stress state is innovatively accounted for in the phase field model by introducing a fictitious strain tensor that is calculated from its corresponding initial stress tensor. It is the sum of the fictitious strain tensor and the strain tensor due to elastic deformation that governs the evolution of the phase field. The simulated results showed that the effect of the swelling pressure (regarded as the initial effective stress for a high swelling clay) on the fracture initiation has been well described by the proposed method. Specifically, the effect of either isotropic or anisotropic stress state on the fracturing process can be well reflected by the phase field approach based on Rankine-type fracture criterion. In contrast, the phase field approach based on the Griffith fracture criterion is more appropriate for the isotropic stress state than the anisotropic stress state because of the Poisson's effect. Moreover, the gas pressure required to trigger the fracturing process needs to exceed the sum of the porewater pressure and the initial stress. The effect of the boundary condition on the evolution of fluid pressure and total stress has been qualitatively captured. It is found that the boundary with higher stiffness leads to a higher gas pressure in the developed fracture and a higher water pressure and total stress in the surrounding porous matrix. In addition, some key experimental findings, such as the preferential gas flow, the build-up of porewater pressure, the almost fully saturated state and the localized consolidation, have been qualitatively captured by the developed phase field model.