Coupled heat transport and fluid flow in porous rocks play a role in many geological phenomena, including the formation of hydrothermal mineral deposits, the productivity of geothermal reservoirs and the reliability of geo-sequestration. Due to the low compressibility of the fluid and rock matrix and the long-time scales the fluid can be treated as incompressible. The solution of the incompressible Darcy flux problem and the advection-dominated heat transport both provide numerically challenging problems typically addressed using methods specialized for the individual equations. In order to avoid the usage of two different meshes and solution approximations for pressure, flux, and temperature we propose to use standard conforming finite-element methods on the same mesh for both problems. The heat transport equation is solved using a linearized finite-element flux corrected transport scheme which introduces minimum artificial diffusion based on the discretized transport problem. The Darcy flux calculation from pressure uses a global post-processing strategy which at the cost of an extra partial differential equation leads to highly accurate flux approximation. In the limit of zero element size the flux is in fact incompressible. We investigate the numerical performance of our proposed method on a test problem using the parallelized modeling environment escript. We also test the approach to simulate convection in geologically relevant scenarios.