The use of characteristics-based methods for the advection-dominated regimes in thermal convection is investigated. An operator-splitting method applied to the advection-diffusion equation for the very large Péclet (Pe) number regime is presented. In this approach two partial differential equations representing both the purely hyperbolic and the parabolic components must be solved simultaneously. This method has been compared with (1) the Galerkin approximation, (2) the streamwise upwinding Petrov-Galerkin method, and (3) the characteristics-based method using the Lagrangian formulation for the time-derivative operator of the advection-diffusion equation. Solution accuracy of the operator-splitting method improves with larger Pe, while the accuracy of other methods deteriorates with Pe. For the nonlinear problem of two-dimensional thermal convection the Lagrangian method is found to be most computationally efficient. With this Lagrangian method, time-dependent, thermal convection solutions of extremely high Rayleigh number (Ra), up to 3 × 10 9, for infinite Prandtl number are obtained. For an aspect ratio of 1.8 the exponent in the scaling of the Nusselt number (Nu) with Ra in time-dependent convection is determined to be 0.301 in the hard turbulent regime and is smaller than in the soft turbulent regime. The behavior of this exponent as a consequence of the transition to hard turbulence agrees with experimental findings. Horizontal Fourier spectra of the thermal fields outside the boundary layers reveal a transition in the high wave-number domain from 1/k to 1/k 2 in the transition from soft to hard turbulent regimes. Analysis of the kinetic energy spectra E(k) shows an asymptotic decay of E(k) close to k -2, for large k, spanning over two decades in wave number for strongly time-dependent convection.