Over the last few decades, numerical models have been proven effective for simulating fracture propagation. Yet very few models are robust enough to simulate real-world fractures. In most of the real-world scenarios, the spatial pattern of fractures is very complicated and two-dimensional models only provide a rough approximation of non-planar fractures. Additionally, in many engineering systems, multiple fractures may co-exist, leading to complex interactions such as crack branching, merging and fragmentation. Quite often, fractures propagate due to the combined effect of thermo-hydro-mechanical loads. As such, modelling real-world fractures remains a formidable challenge. We propose the development of a phase-field based thermo-hydro-mechanical model for fracturing, with a particular emphasis on three-dimensional scenarios. Unlike many existing techniques, a phase-field model inherently incorporates the initiation, propagation, and coalescence of fractures based on energy minimization principles. This feature makes it well-suited for accurately modelling the intricate patterns of cracks in real-world scenarios. Despite this, PFM is mostly limited to 2D scenarios because of its high computational cost, as it requires very fine meshes in the diffused band. We aim to create an efficient phase-field thermo-hydro-mechanical model suitable for analyzing field-scale fracturing in 3D scenarios. The model will be developed in deal.II, a C++ finite element library. The deal.II library allows better performance in computational cost and memory management, making it suitable for large-scale problems. As a first step, we develop a phase-field based fracture simulator to model crack propagation under purely mechanical loading. Benchmark validation examples and preliminary results will be presented in this document, and we'll subsequently expand the model to account for material heterogeneity, applying it to bi-layered and multi-layered formations to showcase the influence of varying material properties on fracture propagation. Thereafter, we shall incorporate MPI and mesh adaptivity into our model to tackle large-scale fracture problems effectively. Finally, we shall develop a fully coupled thermo-hydro-mechanical model to simulate the combined impact of solid deformation, fluid flow and thermal loads.