The advantages of using high-order time integration schemes for thermally coupled flows are assessed numerically. First-, second-, and third-order backward difference schemes are evaluated. The problem is solved in a decoupled manner using a nested iterative algorithm for the Navier–Stokes and energy equations to eliminate decoupling errors. For the space discretization, a stabilized finite element formulation of the variational multiscale type is applied to enable the use of equal order interpolation between the problem unknowns and ensure stable solutions for convection-dominated cases. The integration schemes are compared by solving the flow over a confined square including mixed heat convection in two and three dimensions. Improved numerical approximation of dynamic solutions using high-order schemes is demonstrated in the Richardson number range of 0≤|Ri|≤10 up to a Reynolds number of Re=225.