  This work aims to provide a comprehensive and unified numerical analysis for
non linear system of parabolic variational inequalities (PVIs) subject to
Dirichlet boundary condition. This analysis enables us to establish an
existence of the exact solution to the considered model and to prove the
convergence for the approximate solution and its approximate gradient. Our
results are applicable for several conforming and non conforming numerical
schemes.
