Abstract: A deep understanding of the mechanical properties of gas hydrate-bearing coal (GHBC) is essential for the prevention and control of coal and gas outbursts. The effects of confining pressure, saturation, and heterogeneity on the micro- and meso-mechanical properties of coal are aimed to be evaluated. However, due to the complexity of the internal structure of GHBC, accurately characterizing the quantitative relationship between microscopic heterogeneity and macroscopic mechanical behavior remains a challenge. A DEM based on the Weibull distribution statistical function is proposed in this study to simulate the mechanical properties and failure behavior of GHBC. The validity of the numerical model was validated through laboratory triaxial tests. By varying the heterogeneity within the numerical model, triaxial compression tests were conducted under different confining pressures and saturation levels. The variation patterns of macro- and meso-scale parameters, such as coal strength, contact force chains, and porosity, were analyzed. The experimental results show that: The failure process of both gas-bearing coal (GBC) and GHBC can be divided into the elastic stage, the yield stage, and the strengthening stage. And the stiffness of GHBC is higher than that of GBC, demonstrating superior deformation resistance. As the confining pressures ranges from 12 MPa to 20 MPa, the formation of gas hydrates enhances the bearing capacity of coal specimens. Compared to GBC, the increase rate of peak strength of GHBC is between 9.38% and 36.78%. The numerical simulation results reveal: as the saturation rises from 0 to 80%, the increase rate of peak strength of GHBC is between 8.09% and 67.63%, demonstrating that the presence of hydrates significantly strengthens the load-bearing capacity of coal masses. As the shape factor ranges from 2 to 30, the increase rate of peak strength of GHBC is between 31.62% and 200.94% under different hydrate saturations, and showed a linear increasing trend. As hydrate saturation increases, the influence of the shape factor on specimen strength becomes more pronounced. The proportion of strong force chains increases, the number of cementation failures decreases, and the specimen strength significantly improves. As the hydrate saturation increases from 10% to 80%, the peak value of the average normal contact force fabric tensor rises by 17.03% to 59.27%. Simultaneously, the development of asymmetric “X”-shaped shear bands within the specimen decreases, and the particle rotation velocity gradually reduces. As the shape factor increases, the heterogeneity degree of the specimen weakens, the distribution of strong force chains becomes more uniform, and the structural stability of the specimen system is enhanced. As the shape factor increases from 2 to 30, the peak value of the average normal contact force fabric tensor rises by 7.24% to 31.84%. The rotation speed of particles within the specimen gradually decreases, the width of the shear band gradually decreases, and the degree of damage is mitigated. The axial force is primarily provided by average normal contact forces, while the horizontal force is mainly supported by tangential bonding and friction. The findings facilitate a deeper understanding of its deformation and failure mechanisms, holding significant theoretical importance and engineering application value for enhancing the strengthening effect of hydrates on coal mass.