The Naxos dome, in the middle of the Aegean domain, exposes a part of the former root of the Alpine orogenic belt, and represents a key natural example to investigate the development of gravitational instabilities and their impact on crustal differentiation. The Naxos dome is cored by migmatites with structures depicting second order domes with a diameter of 1 to 2 km nested in the first order deca-kilometer scale dome. Zircon grains from the migmatites have recorded a succession of crystallization-dissolution cycles with a period of 1 to 2 Ma. These features have been attributed to the development of convective and diapiric gravitational instabilities and serve to model the behavior of a former partially molten crust with a volume of fluid method (VOF) that allows large deformation of viscous layers. An originality of this model is that the presence of felsic and mafic domains in a 45 km thick crust of intermediate composition, is mimicked with inclusions of diameter 600 m with contrasted density, viscosity and heat production. Gravitational instabilities are triggered by heating from below and are characterized, first, by the segregation of the inclusions owing to their buoyancy, followed by convection and diapirism. Heavy inclusions accumulate in the lower crust, but are in part entrained in crustal-scale convection cells that control the size of first order domes. Aggregation of the light inclusions at the top of the convective cells leads to the development of second order domes with a diameter of several kilometers that are nested in the first order domes. These models demonstrate the efficiency of gravitational instabilities in the development of domes and in crustal differentiation.