A finite element approach for solving gas-liquid flows with mass transfer and chemical reactions is presented. The incompressible Navier-Stokes equations coupled with a modified k − `` turbulence model for computation of the eddy viscosity are solved by a discrete projection scheme from the family of Pressure Schur Complement methods. Additionally, the simplified two-fluid (Euler-Euler) model based on an analogy of the Boussinesq approximation is adopted for the description of the two-phase flow. A robust solution strategy based on nested iterations is proposed for the numerical treatment of the resulting system of PDEs, which are discretized by high-resolution finite element schemes designed on the basis of algebraic flux correction in order to stabilize the convective terms. The performance of the developed mathematical model is illustrated by a number of numerical results obtained from simulations on three-dimensional unstructured meshes.