This study investigates free thermal convection in heterogeneous porous media and its dependency on grid discretization. Heterogeneous hydraulic conductivity fields are described with a geostatistical approach. Numerical experiments are carried out with the heterogeneous thermal Elder problem. Results are evaluated qualitatively and quantitatively. Key results show that (1) grid convergence for thermal variable-density flow simulations in heterogeneous porous media can be achieved, (2) the degree and structure of a heterogeneous hydraulic conductivity field significantly affect thermal convective flow patterns, (3) steady-state heat distribution depends on the hydraulic conductivity field rather than on thermal conductivity, (4) more variant log hydraulic conductivity fields result in an increase of the heat flux, of the heat stored, and of the upwelling length of thermal plume fingers, (5) an increase in the horizontal,correlation length of the log hydraulic conductivity field leads to an increase of the heat stored and the upwelling length, but it does not affect the heat flux, and (6) an increase in the mean of the log hydraulic conductivity field leads to an increase of the heat flux, of the heat stored, and of the upwelling length. (C) 2016 Elsevier Ltd. All rights reserved.