Abstract: Conventional triaxial compression tests were conducted on gas hydrate-bearing coal(GHBC) with three different saturation levels. Three-dimensional DEM numerical simulations were established based on the anti-rolling model and the parallel bond model respectively, to investigate the influence of saturation and particle shape on the deviatoric stress-strain curves, average coordination number, velocity field, and contact force chains of GHBC. The results indicate that: (1) At saturation levels of 40%, 60%, and 80%, the peak strengths are 14.86, 16.18, and 16.16 MPa respectively, with an increase of 8.7%. The peak strength initially increases with rising saturation but subsequently stabilizes. (2) The anti-rolling contact model demonstrates error rates within 6% for both elastic modulus and peak strength, while the parallel bond model exhibits errors up to 9%. This indicates that both models can characterize the mechanical properties of gas hydrate-bearing coal, although the anti-rolling contact model yields superior accuracy. (3) With increasing saturation, both the initial and final coordination numbers of the specimens increase; the average coordination numbers during loading exhibit distinct patterns between the two sample types; no clear trend emerges in the velocity field evolution; internal contact forces consistently decrease; while peak strength and elastic modulus improve. These changes collectively enhance the mechanical properties of GHBC. (4) Compared to the parallel bond model specimens, the anti-rolling effect enhances effective contacts between particles. This results in force chains that are more numerous and uniformly distributed, primarily forming axial columnar structures, leading to a more stable configuration.