The present work is devoted to the study of unsteady flows of two immiscible viscous fluids separated by a freely moving interface. Our goal is to elaborate a unified strategy for numerical modelling of two-fluid interfacial flows, having in mind possible interface topology changes (like merger or break-up) and realistically wide ranges for physical parameters of the problem. The proposed computational approach essentially relies on three basic components: the finite element method for spatial approximation, the operator-splitting for temporal discretization and the level-set method for interface representation. We show that the finite element implementation of the level-set approach brings some additional benefits as compared to the standard, finite difference level-set realizations. In particular, the use of finite elements permits us to localize the interface precisely, without introducing any artificial parameters like the interface thickness; it also allows us to maintain the second-order accuracy of the interface normal, curvature and mass conservation. The operator-splitting makes it possible to separate all major difficulties of the problem and enables us to implement the equal-order interpolation for the velocity and pressure. Diverse numerical examples including simulations of bubble dynamics, bifurcating jet flow and Rayleigh-Taylor instability are presented to validate the computational method.