In this work, we develop a hydro-damage-mechanical coupled numerical method for hydraulic fracturing in fluid-saturated poroelastic media based on a phase-field regularized cohesive zone model (PF-CZM). By introducing the softening traction-separation laws, the fluid-driven quasi-brittle fracture is simulated and the fracture aperture related to fracture flow is implicitly determined via the crack phase field. The fluid pressure distribution is formulated by a reservoir-fracture fluid field, with the porous reservoir modelled by the Darcy-type flow and the fracture by the modified Darcy-Poiseuille-type flow. The method was first validated by two analytical benchmark problems, including a static Sneddon’s case and a KGD crack propagation problem, and then successfully applied to simulate complex multi-crack fracturing in problems with natural cracks and weakly bonded discontinuities.