Risks and damages associated with lava flows propagation (for instance the most recent Etna eruptions) require a quantitative description of this phenomenon and a reliable forecasting of lava flow paths. Due to the high complexity of these processes, numerical solution of the complete conservation equations for real lava flows is often practically impossible. To overcome the computational difficulties, simplified models are usually adopted, including 1-D models and cellular automata. In this work we propose a simplified 2D model based on the conservation equations for lava thickness and depth-averaged velocities and temperature which result in first order partial differential equations. The proposed approach represents a good compromise between the full 3-D description and the need to decrease the computational time. The method was satisfactorily applied to reproduce some analytical solutions and to simulate a real lava flow event occurred during the 1991--93 Etna eruption.