Abstract. In this work we develop fully discrete (in time and space) numerical schemes for two-dimensional incompressible fluids with mass diffusion, also so-called Kazhikhov-Smagulov models. We propose at most H 1 -conformed finite elements (only globally continuous functions) to approximate all unknowns (velocity, pressure and density), although the limit density (solution of continuous problem) will have H 2 regularity. A backward Euler in time scheme is considered decoupling the computation of the density from the velocity and pressure.Unconditional stability of the schemes and convergence towards the (unique) global in time weak solution of the models is proved. Since a discrete maximum principle cannot be ensured, we must use a different interpolation inequality to obtain the strong estimates for the discrete density, from the used one in the continuous case. This inequality is a discrete version of the GagliardoNirenberg interpolation inequality in 2D domains. Moreover, the discrete density is truncated in some adequate terms of the velocity-pressure problem.