Numerical experiments of the development of gravitational instabilities in 3D view (top) and 2D view (bottom). The ambient medium contains layers that are less dense and less viscous (white, felsic material) or denser and more viscous (black, mafic material). Basal heating is applied, isotherms are depicted by labeled colored solid lines and viscosity is indicated in the blue to red color bar. Viscosity of the modeled crust varies from 10¹⁶ to 10²³ Pa.s as a function of temperature and strain rate, while density decreases with temperature. This setting mimics a partially melting crust with felsic and mafic layers. After 2.0 Myr, the evolution of the model during heating is characterized by redistribution of material phases according to their buoyancy, namely upward motion of the white material and sinking of the black material. Composition-driven diapirs of buoyant white material develop at ca. 2.5 Myr and ca. 5 Myr, after their accumulation below the dense and more viscous black layers. Thermal-convection starts after 10 Myr and the size of convection cells increases with the thickness of low-viscosity crust (less than ∼ 10¹⁹ Pa.s), to reach about 25 km in diameter at ca. 20 Myr. Some of the black and white materials are entrained in convection but most of the dense material accumulates at the bottom of the modeled crust, whereas the buoyant material continues to accumulate at the top of the convection cells, forming diapirs with a ∼ 5 km diameter. An animation of this simulation is available here: https://www.youtube.com/playlist?list=PLCdeVmWEruk1x-kyWxAkl4be7b2yw5eG9