Rigorous coupling of fracture–porous medium fluid flow and topologically complex fracture propagation is of great scientific interest in geotechnical and biomechanical applications. In this paper, we derive a unified fracture–porous medium hydraulic fracturing model, leveraging the inherent ability of the variational phase-field approach to fracture to handle multiple cracks interacting and evolving along complex yet, critically, unspecified paths. The fundamental principle driving the crack evolution is an energetic criterion derived from Griffith’s theory. The originality of this approach is that the crack path itself is derived from energy minimization instead of additional branching criterion. The numerical implementation is based on a regularization approach similar to a phase-field model, where the cracks location is represented by a smooth function defined on a fixed mesh. The derived model shows how the smooth fracture field can be used to model fluid flow in a fractured porous medium. We verify the proposed approach in a simple idealized scenario where closed form solutions exist in the literature. We then demonstrate the new method’s capabilities in more realistic situations where multiple fractures turn, interact, and in some cases, merge with other fractures.