<p>This study presents a phase-field modeling framework that combines an interfacial phase-field approach with adaptive mesh refinement to simulate thermomechanically-induced fractures in layered rocks. Meanwhile, this study implements the interfacial phase-field method in COMSOL. The interfacial phase-field method captures smooth transitions in material properties across bedding planes, avoiding explicit interface modeling while accurately representing mechanical and thermal responses near these interfaces. To reduce the computational cost typical of phase-field fracture simulations, an adaptive mesh refinement strategy is employed using the COMSOL API. The mesh refinement is dynamically guided by the phase-field variable within COMSOL’s Application Builder, enabling focused refinement around evolving cracks while maintaining coarser meshes elsewhere. The coupled four-field system (temperature, displacement, phase-field, and interfacial phase-field) is solved through segregated solution steps (staggered solution scheme). Compared with uniform mesh refinement, the adaptive approach significantly reduces computational demands without sacrificing accuracy in predicting crack paths and fracture morphologies. Validation through multiple numerical examples under quasi-static thermomechanical loading demonstrates the framework’s capability to capture complex fracture processes influenced by thermal effects and bedding-plane heterogeneity. This approach offers a robust and efficient tool for modeling fractures in layered rocks, with practical implications for geothermal energy extraction, nuclear waste disposal, and deep underground engineering.</p>