We introduce a time discretization for the incompressible Navier-Stokes equations with variable density and thermal convection. These equations governs the dynamics of many problems in geosciences and engineering such as ocean-atmospheric interactions, convection in the Earth's mantle layers, production of alumina with Hall-H\'{e}roult cell or storage of energy using liquid metal batteries where the Maxwell equations can also be involved for some of the above problems. The main difficulties of such setups consists of tracking the interface between fluids (with different density, viscosity, thermal conductivity, etc.) and approximating the Navier-Stokes equations whose system involves a diffusive operator and a time derivative operator that both have space-time dependent coefficients. We note that similar issues arise with the temperature equation where the heat capacity and thermal diffusivity also are space-time dependent variable. To tackle such problems, we propose a numerical method that uses the momentum and internal energy as primary variables. The interface between the fluids is tracked using a level set technique. The incompressibility constraints is enforced using an artificial compressibility technique which, unlike projection method, does not require to solve an elliptic equation to update the pressure. Thanks to an appropriate treatment of the diffusive operators and nonlinear terms, the resulting algorithm only uses time-independent stiffness matrices that can be assembled and preconditioned at initialization. Moreover, the proposed method is suitable for high order finite element and spectral discretization. We establish the stability and time-convergence properties of the proposed first-order algorithm. Applications to the generation of turbulent thermal convection in two layered systems in three dimensional domains will be presented.