The numerical modeling of compressible multiphase flows is of high interest for several engineering applications. In this work, we focus on the study of pyroclastic flows arising from volcanic eruptive events. An accurate evaluation of the effects of this multiphase flow is of crucial importance for the Civil Protection for the preservation of urban settlements in volcanic areas. In this work, we propose a Finite Element formulation for the simulation of pyroclastic flows at conditions of thermal and kinetic equilibrium. This analysis belongs to the class of advection dominated problems, which are known to suffer from numerical instabilities. The required stabilization is provided by using a Variational Multiscale Method, typically used for monophase compressible flows and extended here for the first time to multiphase flows. The stabilized formulation is validated with benchmark problems for compressible flows, which are solved for both monophase and multiphase cases. A throughout comparison of the numerical results of the two different flows is also presented. Moreover, the numerical formulation is applied to the simulation of representative cases of pyroclastic flows considering large-scale computational domains and realistic material properties and initial thermal-kinematic conditions. The numerical analyses presented show the accuracy of the proposed method for the simulation of compressible multiphase flows and its suitability for risk assessment studies of urban settlements prone to be affected by pyroclastic gravity currents.