<p>To deepen understanding of the hydraulic fractures (HFs) initiation and propagation in coal with complex discrete fracture networks, a fully coupled hydro-mechanical-damage (HMD) continuum-discrete model was developed, incorporating the dual-medium characteristics of coal. The model represents the matrix system as a continuum medium system and the fracture system as a discrete medium system. Field-scale simulations of hydraulic fracturing in naturally fractured coal demonstrated the model's ability to capture direct HM coupling effectively, indirect coupling processes involving property alterations in the matrix and natural fractures (NFs), and interactions between the HFs and NFs. The results reveal that the cumulative damage area in the matrix increases almost linearly during the fracturing process, while the fluid pressure and effective mean stress in the HFs exhibit rapid growth, sharp decline, and eventual stabilization. On both sides of the HFs, the vertical displacement of the matrix shows discontinuous variations. The fracture system, as the primary flow pathway, dominates the fluid pressure distribution within the matrix and significantly influences the propagation of the HFs. The evolution of the ratio of shear-dilating NF length to the total NF length exhibits a trend of initial increase followed by a decrease under HMD coupling conditions. Notably, the HMD coupling effect significantly enhances the shear dilation of NFs compared with the HM coupling effect. This model provides substantial advantages in capturing the complex HMD coupling effects within the dual-medium system, offering novel insights into the evolution of stress, fluid pressure, and fracture propagation during hydraulic fracturing in naturally fractured coal.</p>