A numerical model of two-dimensional Rayleigh-B¿nard convection is used to study the relationship between the surface heat flow (or Nusselt number) and the viscosity at the base of the lithosphere. Newtonian or non-Newtonian, temperature- and pressure-dependent rheologies are considered. In the high Rayleigh number time-dependent regime, calculations yield Nu∝RaBL1/3beff-4/3, where beff is the effective dependence of viscosity with temperature at the base of the upper thermal boundary layer and RaBL is the Rayleigh number calculated with the viscosity &ngr;BL (or the effective viscosity) at the base of the upper thermal boundary layer. The heat flow is the same for Newtonian and non-Newtonian rheologies if the activation energy in the non-Newtonian case is twice the activation energy in the Newtonian case. In this chaotic regime the heat transfer appears to be controlled by secondary instabilities developing in thermal boundary layers. These thermals are advected along the large-scale flow. The above relationship is not valid at low heat flow where a stationary regime prevails and for simulations forced into steady state. In these cases the Nusselt number follows a trend Nu∝RaBL1/5beff-1 for a Newtonian rheology, as predicted by the boundary layer theory. We argue that the equilibrium lithospheric thickness beneath old oceans or continents is controlled by the development of thermals detaching from the thermal boundary layers. Assuming this, we can estimate the viscosity at the base of the stable oceanic lithosphere. If the contribution of secondary convection to the surface heat flux amounts to 40 to 50 mW m-2, the asthenospheric viscosity is predicted to be between 1018 and 2¿1019 Pa s. ¿ 1999 American Geophysical Union