Interfaces, free or internal, greatly influence the physical properties and stability of materials microstructures. Of particular interest are the processes that occur due to anisotropic interfacial properties. In the case of grain boundaries (GBs) in metals, several experimental observations revealed that an initially flat GB may facet into hill-and-valley structures with well defined planes and corners/edges connecting them. Herein, we present a diffuse interface model that is capable of accounting for strongly anisotropic GB properties and capturing the formation of hill-and-valley morphologies. The hallmark of our approach is the ability to independently examine the various factors affecting GB faceting and subsequent facet coarsening. More specifically, our formulation incorporates higher order expansions to account for the excess energy due to facet junctions and their non-local interactions. As a demonstration of the modeling capability, we consider the Σ5 〈001〉 tilt GB in body-centered-cubic iron, where faceting along the {210} and {310} planes was experimentally observed. Atomistic calculations were utilized to determine the inclination-dependent GB energy, which was then used as an input in our model. Linear stability analysis and simulation results highlight the role of junction energy and associated non-local interactions on the resulting facet length scales. Broadly speaking, our modeling approach provides a general framework to examine the microstructural stability of polycrystalline systems with highly anisotropic GBs.