We develop an algorithm to solve numerically the Navier–Stokes equations using a finite element method. The algorithm uses a fractional step approach that allows the use of equal interpolation spaces for the pressure and velocity fields and can be used to solve both compressible and incompressible flows. The standard Galerkin method is used to space-discretize the equations in which convective terms are dominant because the equations are reformulated in a characteristics co-moving frame providing thus the required artificial diffusion in a consistent way. The algorithm depends on four different parameters. Depending on their values, a fully implicit, a semi-implicit or a fully explicit solution can be obtained. The proposed algorithm may be useful in solving a wide spectrum of problems in geology, engineering and geosciences. Several flows frequently encountered in practical applications such as incompressible, slightly compressible or perfect gas are taken into account. As an example, we use it to model the withdrawal of magma from shallow chambers during explosive volcanic eruptions. Our model constitutes a first attempt to characterize the temporal evolution of the most relevant physical parameters during such a process.