The phase-field method has emerged as a sophisticated technique for simulating crack initiation, propagation, and coalescence. This approach employs a damage field, termed the phase field, to represent the material's state from intact to fully fractured. The phase-field approach is well known to yield a highly nonlinear and non-convex system of equations. Therefore, the design of efficient and robust solution methods to address such challenging systems of equations is of critical importance. In this poster, we present novel scalable algorithms, preconditioning strategies, and high-performance implementation of a finite-element-based solver for the phase-field fracture formulation. We provide novel insights into fracture-infilling mechanisms of sedimentary layers and illustrate a geological benchmark for the phase-field community.