The efficiency of implicit and semi-implicit time integration codes based on backward differentiation and extrapolation formulas for the solution of the three-dimensional Boussinesq thermal convection equations in rotating spherical shells was studied in [5] at weakly supercritical Rayleigh numbers R, moderate (10−3) and low (10−4) Ekman numbers, E, and Prandtl number = 1. The results presented here extend the previous study and focus on the effect of and R by analyzing the efficiency of the methods for obtaining solutions at E = 10−4, = 0.1 and low and high supercritical R. In the first case (quasiperiodic solutions) the decrease of one order of magnitude does not change the results significantly. In the second case (spatio-temporal chaotic solutions) the differences in the behavior of the semi-implicit codes due to the different treatment of the Coriolis term disappear because the integration is dominated by the nonlinear terms. As in [5], high order methods, either with or without time step and order control, increase the efficiency of the time integrators and allow to obtain more accurate solutions.