The origin of stable old continental cratonic roots is still debated. We present numerical modelling results which show rapid initial formation during the Archaean of continental roots of ca. 200 km thick. These results have been obtained from an upper mantle thermal convection model including differentiation by pressure release partial melting of mantle peridotite. The upper mantle model includes time-dependent radiogenic heat production and thermal coupling with a heat reservoir representing the Earth's lower mantle and core. This allows for model experiments including secular cooling on a time-scale comparable to the age of the Earth. The model results show an initial phase of rapid continental root growth of ca. 0.1 billion year, followed by a more gradual increase of continental volume by addition of depleted material produced through hot diapiric, convective upwellings which penetrate the continental root from below. Within ca. 0.6 Ga after the start of the experiment, secular cooling of the mantle brings the average geotherm below the peridotite solidus thereby switching off further continental growth. At this time the thickness of the continental root has grown to ca. 200 km. After 1 Ga of secular cooling small scale thermal instabilities develop at the bottom of the continental root causing continental delamination without breaking up the large scale layering. This delaminated material remixes with the deeper layers. Two more periods, each with a duration of ca. 0.5 Ga and separated by quiescent periods were observed when melting and continental growth was reactivated. Melting ends at 3 Ga. Thereafter secular cooling proceeds and the compositionally buoyant continental root is stabilized further through the increase in mechanical strength induced by the increase of the temperature dependent mantle viscosity. Fluctuating convective velocity amplitudes decrease to below 10 mma−1 and the volume average temperature of the sub-continental convecting mantle has decreased ca. 340 K after 4 Ga. Surface heatflow values decrease from 120 to 40 mW m−2 during the 4 Ga model evolution. The surface heatflow contribution from an almost constant secular cooling rate was estimated to be 6 mW m−2, in line with recent observational evidence. The modelling results show that the combined effects of compositional buoyancy and strong temperature dependent rheology result in continents which overall remain stable for a duration longer than the age of the Earth. Tracer particles have been used for studying the patterns of mantle differentiation in greater detail. The observed ( p, T, F, t)-paths are consistent with proposed stratification and thermo-mechanical history of the depleted continental root, which have been inferred from mantle xenoliths and other upper mantle samples. In addition, the particle tracers have been used to derive the thermal age of the modelled continental root, defined by a hypothetical closing temperature.