Crack nucleation in thermal-barrier coating (TBC) systems subjected to a monotonic cooling process is studied. The TBC system is modeled using the finite element method, where cracks are represented as discrete discontinuities across continuum elements using the partition-of-unity method. The numerical implementation used for crack nucleation is based on an algorithm where, at insertion of a discontinuity, the traction response is derived from a cohesive zone model that has been modified to (i) behave like an initially rigid cohesive model, and to (ii) ensure smoothness of the traction-separation law at zero crack opening. Accordingly, an adequate convergence behavior of the numerical formulation can be warranted in boundary value problems of systems with relatively complex geometries. In the present numerical study, a comparison is made between TBC systems composed of different constitutive models. The fracture patterns and evolutions of the overall crack growth resulting from the simulations clearly illustrate the importance of accounting for the effects of plasticity in the bond coating and anisotropy in the top coating. The computed fracture profile is in good correspondence with experimental observations reported in the literature.