The Naxos dome, in the middle of the Aegean domain, exposes the former root of the Alpine orogenic belt and represents a key natural example to investigate the development of gravitational instabilities during orogenic evolution and their impact on crustal differentiation. The Naxos dome is cored by migmatites with structures depicting second order domes with a diameter of 1–2 km nested in the first order deca-kilometer scale dome that formed at the onset of orogenic collapse. Zircon grains from the migmatites record a succession of crystallization-dissolution cycles with a period of 1–2 Myr. These features have been attributed to the development of convective and diapiric gravitational instabilities, related to thermally induced and compositional buoyancy. In this paper, we test the pertinence of this model with a thermal-mechanical numerical experiment performed with a volume of fluid method (VOF) known to preserve material phase interfaces during large deformation of viscous layers. Partial melting of the crust is modeled by strain-rate and temperature dependent viscosity and temperature dependent density. Moreover, horizontal layers with density, viscosity and heat production variations mimic more felsic or more mafic lithologies in a crust of intermediate composition. With basal heating, gravitational instabilities initiate with local segregation of the buoyant versus heavier layers, followed by diapiric upwelling of buoyant pockets of aggregated less dense material. Convection starts after 5 Myr, approximately when half of the crust has a viscosity lower than 1019 Pa s. The size of the convection cells increases as the temperature rises in the crust and reaches ∼25 km in diameter after ca. 20 Myr, which defines the size of first order domes. Some of the heterogeneous material is entrained in the convection cells with a revolution period of 1 to 3 Myr. However, most of the denser material accumulates in the lower crust, while the buoyant material segregates at the top of the convection cells and forms diapirs that correspond to second order domes, of several kilometers in diameter and nested within the first order domes. This model, which reproduces the first order characteristic dimensions of the Naxos nested domes and the periodicity of their zircon geochronological record, demonstrates the efficiency of gravitational instabilities in the formation of migmatite domes and, more generally, in the multi-scale dynamics of crustal differentiation leading to a felsic upper crust, an intermediate middle crust and a mafic lower crust.