Development of a computational solver based upon the high-order, high-resolution Galerkin finite element method (FEM) often encounters two challenges: First, the Galerkin FEMs give central approximations to the differential operators and their use in the simulation of the convection-dominated flows may lead to the dispersion errors yielding entirely wrong numerical solutions. Secondly, high-order, high-resolution numerical methods are known to produce high wave-number oscillations in the vicinity of shocks/discontinuities in the numerical solution adversely affecting the stability of the method. We present the stabilized finite element method for plasma fluid models to address the two challenges. The numerical stabilization is based on two strategies: Variational Multiscale (VMS) and the shock-capturing approach. The former strategy takes into account (the approximation of) the effect of the unresolved scales onto resolved scales to introduce upwinding in the Galerkin FEM. The latter adaptively adds the artificial viscosity only in the vicinity of shocks. These two strategies can be used to improve the stability and robustness of the numerical methods to solve a wide range of physical problems in fluid dynamics, plasma physics, astrophysics, etc. In this work, we demonstrate the use of the stabilized FEM to perform simulations of complex plasma dynamics in tokamaks. Massive material injection (MMI) experiments in the tokamak plasma constitute gas-plasma interactions (GPI) that can be described using magnetohydrodynamics (MHD) equations along with the impurities transport. The physical phenomena involved are convection-dominated, anisotropic and highly nonlinear containing shocks and strong local sources. The numerical stabilization strategies developed in this work are applied to the high-order bi-cubic Hermite Bézier FEM in the computational framework of the nonlinear MHD code JOREK. The stabilized FEM is then used to simulate the complex MHD dynamics occurring during GPI in tokamak plasma. The physical models and numerical methods presented in this work will be used to perform large scale MHD simulations to help understand plasma instabilities that are critical for the design of tokamak elements.