Abstract
This work aims to provide a comprehensive and unified numerical analysis for a nonlinear system of parabolic variational inequalities (PVIs) subject to Dirichlet boundary condition. This analysis enables us to establish the existence of an 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 nonconforming numerical schemes.