The aim of the paper is twofold. On one hand the authors want to present a new technique called \(p\)-caloric approximation, which is a proper generalization of the classical compactness methods first developed by DeGiorgi with his Harmonic Approximation Lemma. This last result, initially introduced in the setting of Geometric Measure Theory to prove the regularity of minimal surfaces, is nowadays a classical tool to prove linearization and regularity results for vectorial problems. Here the authors develop a very far reaching version of this general principle devised to linearize general degenerate parabolic systems. The use of this result in turn allows the authors to achieve the subsequent and main aim of the paper, that is, the implementation of a partial regularity theory for parabolic systems with degenerate diffusion of the type \(\partial_t u - \mathrm{div} a(Du)=0\), without necessarily assuming a quasi-diagonal structure, i.e. a structure prescribing that the gradient non-linearities depend only on the the explicit scalar quantity.