Polygonal faults are thought to potentially compromise the integrity of regional caprocks and are widely developed in mudrocks that are targeted for unconventionals exploration. Polygonal Fault Systems (PFS) contain networks of exclusively normal faults that intersect bedding planes at a wide variety of azimuths. This fact alone suggests a non-tectonic origin and recent arguments focus on a constitutive control on PFS formation. Geomechanical forward modelling has the potential to shed light on the evolution of material and stress state as geological structures form. The evolution of polygonal faults is studied here using the finite strain forward modelling code ELFEN FM. More specifically the recently suggested diagenetic trigger for PFS formation is investigated, and the process of incorporating this into a critical-state based model is outlined and discussed. The results demonstrate the formation of PFS in both 2D plane-strain and full 3D simulations as validation of both the genetic argument and the computational approach. PFS are additionally observed to be sensitive to subtle horizontal stress anisotropy e.g. around salt diapirs or larger tectonic faults. As such, there is the potential for these fault networks to act as 'paleo-stress piezometers'. This study investigates this by incorporating subtle stress anisotropy into the model and observing the influence on fault alignment. Results indicate that even for subtle stress anisotropy, preferential fault alignment is observed and this is consistent with observation.