# Formation of the Naxos nested domes and crustal differentiation by convection and diapirism

Aurélie Louis-Napoléon<fnref n="1" />,<fnref n="2" />, Olivier Vanderhaeghe<fnref n="2" />,*, Muriel Gerbault<fnref n="2" />, Roland Martin<fnref n="2" /> and Thomas Bonometti<fnref n="1" />

Received: 19 April 2023 / Accepted: 12 August 2024 / Publishing online: 15 October 2024

**Abstract** – 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 10¹⁹ 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.

**Keywords:** Geology of Naxos / Aegean domain / Nested domes / Gravitational instabilities / Crustal convection / Diapirism

**Résumé** – **Formation des dômes emboités de Naxos par convection et diapirisme.** Le dôme de Naxos, au milieu du domaine égéen, expose la racine de la ceinture orogénique alpine et représente une cible naturelle de 1er choix pour aborder le développement d'instabilités gravitaires au cours de l'évolution orogénique et leur impact sur la différenciation crustale. Le cœur du dôme de Naxos est constitué de migmatites dont la structure souligne des dômes de second ordre avec un diamètre de l'ordre de 1 à 2 km emboités dans le dôme de 1er ordre de taille déca kilométrique qui se sont formés au début de l'effondrement orogénique. Les grains de zircon des migmatites ont enregistré une succession de cycles de cristallisation-dissolution avec une périodicité de l'ordre de 1 à 2 Ma. Ces caractéristiques sont attribuées au développement d'instabilités gravitaires convectives et diapiriques induites par des différences de densité

<footnote n="1">¹ Institut de Mécanique des Fluides de Toulouse (IMFT), Université de Toulouse, CNRS, Toulouse INP, UPS, Toulouse, France</footnote>

<footnote n="2">² Géosciences Environnement Toulouse (GET), Université de Toulouse, CNRS, UPS, IRD, CNES, Toulouse, France</footnote>

<footnote>*e-mail: olivier.vanderhaeghe@get.omp.eu</footnote>
liées à la température et à la composition. Dans ce papier, nous testons la pertinence de ce modèle avec une expérience numérique thermomécanique conduite avec la méthode « volume of fluid » (VOF) qui assure un suivi des interfaces entre phases matérielles au cours de la déformation de milieux visqueux. La fusion partielle est modélisée par une viscosité dépendante de la vitesse de déformation et de la température et par une densité fonction de la température. La présence de niveaux horizontaux présentant des contrastes de densité, de viscosité et de production de chaleur par rapport au milieu ambiant mimique le comportement d'une croûte de composition intermédiaire comprenant des niveaux felsiques et mafiques. Le réchauffement imposé à la base du modèle induit le développement d'instabilités gravitaires qui débutent par la ségrégation locale des matériaux en fonction de leur densité puis se poursuit par l'accumulation de proches de matériel peu dense et la formation de diapirs. La convection débute après 5 Myr, alors qu'environ la moitié de la croûte modélisée a une viscosité plus faible que $10^{19}$ Pa.s. La taille des cellules de convection augmente avec la température et atteint un diamètre de ~25 km après ca. 20 Myr, ce qui définit la taille des dômes de 1er ordre. Une partie du matériel est entraîné dans la convection, indépendamment de sa densité avec une période de révolution de 1 à 3 Myr. Cependant, la majorité du matériel dense est accumulée à la base de la croûte alors que le matériel moins dense est ségrégé au sommet des cellules de convection et forme des diapirs qui correspondent à des dômes de 2nd ordre avec un diamètre de quelques kilomètres, enveloppés dans les dômes de 1er ordre. Ce modèle, qui reproduit la taille caractéristique des dômes emboîtés de Naxos et la périodicité de l'enregistrement géochronologique des grains de zircon qu'ils contiennent, démontre l'efficacité des instabilités gravitaires pour la formation des dômes migmatitiques et plus généralement de la dynamique multi-échelle de la différenciation crustale conduisant à la formation d'une croûte litée avec une croûte supérieure felsique, une croûte moyenne de composition intermédiaire et une croûte inférieure mafique.

**Mots-clés** : Géologie de Naxos / Domaine Egéen / Dômes emboîtés / Instabilités gravitaires / Convection crustale / Diapirisme

# 1 Introduction

The formation of domes cored by migmatites is a long standing debated topic since their first description in the Baltic Shield (Eskola, 1948). Some authors have insisted on the role of buoyancy in the development of gravitational instabilities (Brun et al., 1981; Collins, 1989; Ramberg, 1981; Talbot, 1979) while others have invoked the role of tectonic forces in compressional contexts (Burg and Podladchikov, 2000; Myers and Watkins, 1985; Porada and Berhorst, 2000) or in extensional context related to metamorphic core complexes (Brun et al., 1994; Buck, 1991; Coney and Harms, 1984; Davis, 1983; Le Pourhiet et al., 2012). These different scenarios are not mutually exclusive but might be distinguished on the basis of structural, metamorphic and geochronological data (Brun, 1983; Burg et al., 2004; Van Kranendonk et al., 2004; Whitney et al., 2004; Yin, 2004). On the other hand, melt/solid segregation has also been invoked as a major mechanism to redistribute chemical elements, leading to crustal differentiation (Brown, 2001; Cruden et al., 1995; Petford et al., 2000; Sawyer, 1994; Vanderhaeghe, 2001; Vanderhaeghe et al., 2009).

Naxos Island, in the central part of the Aegean domain, provides exceptional exposure of a dome cored by migmatites and has been an emblematic target to discuss the driving forces for dome formation. The first geological map of Naxos (Jansen, 1973) formed the basis for interpreting the migmatite-cored dome, and associated concentric metamorphic isograds, as a diapir (Jansen and Schuiling, 1976). The identification of a major detachment juxtaposing a migmatite-bearing footwall of high-grade metamorphic rocks, to a hanging-wall comprising low-grade metamorphic rocks associated with syntectonic sedimentary deposits, led to a revised interpretation of dome formation as a consequence of crustal extension (Buick, 1991; Gautier et al., 1993; Lister et al., 1984; Rey et al., 2011). Other authors have proposed that the Naxos dome formed due to fold interferences during N-S extension associated with E-W shortening

(Avigad et al., 2001; Buick, 1991; Lamont et al., 2020, 2023; Linnros et al., 2019; Peillod et al., 2017, 2021a) or as a part of strike slip structure (Le Pourhiet et al., 2012). Detailed structural analysis of the migmatitic rocks of Naxos revealed the presence of second order domes nested within the main dome, which have been attributed to gravitational instabilities (Kruckenberg et al., 2011; Vanderhaeghe et al., 2018). This interpretation is consistent with thermo-mechanical models of either diapirism or convection of the partially molten crustal root when it is heated from below. It was applied to the Aegean domain (Schenker et al., 2012) but also to other natural examples (Babeyko et al., 2002; Cruden et al., 1995; Riel et al., 2016; Schmeling et al., 2019; Weinberg, 1997; Weinberg and Schmeling, 1992; Zuza and Cao, 2023). Nevertheless, the development of gravitational instabilities does not preclude the implication of tectonic forces. Indeed, the geological record of Naxos might reflect a combination of gravitational collapse, lateral extrusion and gravitational instabilities along a retreating convergent plate boundary (Gautier et al., 1999; Jolivet and Brun, 2010, 2021; Vanderhaeghe and Teyssier, 2001; Vanderhaeghe et al., 2007).

In this contribution, we further explore the role of gravitational instabilities (thermally-driven convection and compositionally-driven diapirism) on the formation of nested migmatitic domes based on the example of Naxos. In previous work, we have first investigated analytically the conditions for which a continental felsic crust may convect, when it is either heated from below or subject to internal heating; dimensional analysis provided realistic critical crustal thickness and average viscosities (Vanderhaeghe et al., 2018, cf. Appendix Fig. A1). Then we investigated with numerical models the behavior of a low-viscosity crust resulting from partial melting. We employed a Volume of Fluid method (VOF), which allows us to efficiently track material interfaces in contexts of very large deformation (Louis-Napoléon et al., 2020, 2022). In these papers we tested the influence of several viscosity evolution laws as a function of temperature and
several melt fraction thresholds. We also explored the impact of heterogeneities, termed “inclusions”, of different sizes and of contrasted viscosity and density relative to the ambient medium. These simulations allowed to identify four different flow regimes that describe the motion of the inclusions : (i) a “segregation” regime, (ii) a “diapiric” regime controlled by the buoyancy of the inclusions relative to the ambient medium, (iii) a “suspension” regime whereby thermally-driven convection drags the inclusions within crustal scale convective cells, and (iv) a “layering” regime whereby convection is concomitant with the accumulation of inclusions at the bottom and at the top of the convective cells for the denser and lightest ones, respectively. In this regime, compositionally-driven diapirs tend to form second order domes at the top of the convective cells (Appendix Fig. A2). For a full description of the methods, benchmarks and parametric investigation, the reader is referred to Louis-Napoléon et al. (2020, 2022, 2024).

In the present paper we present the result of a thermal-mechanical model designed to investigate the effect of gravitational instabilities in a partially molten crust, neglecting the role of tectonic forces. This is clearly an oversimplification and discrepancies between the model and Naxos geological record allow us to further discuss the impact of the regional geodynamic context and its influence on dome formation. We pre-selected the physical parameters (thermal evolution, viscosity and density contrasts of heterogeneities) required to reproduce the characteristic length and time scales for the formation of the nested migmatite domes of Naxos Island, determined from the structural, petrological and geochronological data available for Naxos. The goal here is not to reproduce the parametric analysis performed in Louis-Napoleon et al. (2020, 2022, 2024), but rather to expose the mechanical feasibility for the congruent development of both convection and diapiric structures, leading to the formation of nested domes by gravitational instabilities in the specific case of Naxos.

## 2 Naxos Island geological context in the frame of the Hellenides-Aegean

The Hellenides-Aegean orogenic belt formed due to convergence between Africa and Eurasia starting in the Middle Cretaceous (Dewey and Sengör, 1979; Dercourt et al., 1986) and it has been marked by slab retreat since at least Miocene times (Fytikas et al., 1984; Spakman et al., 1988). Tectonic-petrologic reconstructions lead to at least 45 km thick orogenic crust that would have developed prior to gravitational collapse (Gautier et al., 1999; Jolivet and Brun, 2010; Ring et al., 2010; Vanderhaeghe and Teyssier, 2001) compared to the present-day crustal thickness of 20–25 km (Tirel et al., 2004).

The most complete geological section of the Aegean domain is exposed on Naxos (Gautier et al., 1993; Lamont et al., 2020; Peillod et al., 2021b; Vanderhaeghe et al., 2007). The major rock units in the island represent a metamorphic core complex with (i) an Upper Unit of low-grade metamorphic rocks including an ophiolitic mélange and Miocene coarse-grained sediments, (ii) a Middle Unit dominated by a series of micaschist and marble layers including amphibolite boudins, attributed to a Mesozoic continental margin sedimentary sequence, which is part of the Cycladic Blueschist

Unit, and (iii) a Lower Unit made of migmatites (Fig. 1). The contact between the Upper Unit and Middle Unit is marked by a low-angle detachment bearing a dominant NNE-SSW trending stretching lineation (Buick, 1991; Gautier et al., 1993; Urai et al., 1990). The timing of the detachment fault activity is bracketed between 12 and 9 Ma by (i) emplacement of a syntectonic granodiorite pluton that is intrusive in the Middle Unit and marked by a mylonitic to cataclastic fabric with kinematic criteria consistent with a top to the NNE sense of shear (Brichau et al., 2006; Gautier et al., 1993; Keay et al., 2001; Seward et al., 2009) and (ii) deposition of syn to post-tectonic Miocene clastic sediments that are cross-cut by normal faults and locally display a fan shape geometry (Gautier et al., 1993; Vanderhaeghe et al., 2007). The Middle Unit is characterized by the superposition of F₁ and F₂ isoclinal folds that result in the development of a composite foliation bearing a NNE-SSW trending stretching lineation (Buick, 1991). Mineral relics of blueschist facies HP/LT metamorphism of the Middle Unit have been identified along the southern tip of the Island (Avigad, 1998) and dated at ca. 50 Ma by argon thermochronology on white mica (Wijbrans and McDougall, 1988). These rocks are overprinted by a MP/MT metamorphic event grading from greenschist to amphibolite facies dated between 35 and 20 Ma (Buick and Holland, 1989; Duchêne et al., 2006; Jansen and Schuiling, 1976; Keay et al., 2001; Martin et al., 2006, 2008; Peillod et al., 2017).

The contact between the Middle and Lower Unit corresponds to the melt-in isograd, which is gradual and is interpreted as a MP/MT metamorphic gradient superimposed on the former HP/LT gradient. It is further marked by tilting and transposition of the metamorphic foliation into a syn-migmatic foliation marked by the alternation of leucosome and mesosome +/− melanosome, which delineates the core of the Naxos dome (Kruckenberg et al., 2011; Vanderhaeghe, 2004). The migmatites are dominated by diatexites made of a heterogeneous granite matrix enclosing numerous enclaves of metapelite, marble and amphibolite. Diatexites are coring the nested domes and are locally surrounded by layers of metatexitic paragneiss and/or marble that can reach a thickness of several hundred meters. Overall, diatexites are dominated by quartz and feldspar and have a granodioritic composition while the metamorphic series of the Middle Unit are dominated by micaschists composed of micas and quartz alternating with calcite marble and amphibolite. In other words, the diatexites are on average more felsic than the metamorphic series and this might reflect differentiation of the partially molten crust.

The Naxos dome has an elliptical shape with a long axis of ~12 km and a short axis of ~5 km. The migmatites fabric has been identified both in the field and by measurements of the magnetic susceptibility of oriented rock samples: they define concentric foliation trajectories locally associated with a radial lineation that delimit second order domes with diameters of 2–3 km (Kruckenberg et al., 2010, 2011; Vanderhaeghe, 2004; Fig. 1).

The Naxos dome is also delineated by concentric metamorphic isograds. Kyanite-bearing micaschist at the contact between the Middle Unit and Lower Unit, located structurally above the melt+ isograd, have recorded a pressure of ~0.7 GPa and a temperature of ~650 °C, which are consistent with the location of the granitic solidus at about 20 km depth. If we add the present day crustal thickness of
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Geology of the Naxos metamorphic core complex characterized by nested domes within a larger elliptical dome cored by migmatites (modified after Kruckenberg *et al.*, 2011).</figcaption></figure>

25 km, this leads to a thickness of about 45 km at the time of the metamorphic pressure peak (Duchêne *et al.*, 2006). On the other hand, the sillimanite-bearing migmatites of the Lower Unit have recorded a metamorphic pressure of ~ 0.4 GPa, which is 0.3 GPa less than the pressure recorded by the nearby overlying kyanite-bearing micaschist. This apparent paradox is resolved if we consider that 9 km of unroofing by crustal extension occurred between the time at which the maximum crustal thickness was reached, recorded at 0.7 GPa in the kyanite-bearing micaschist, and the time of final crystallization of the partially molten crust recorded in the sillimanite-bearing migmatites (Duchêne *et al.*, 2006).

Zircon U-Pb ages obtained from the migmatites of the Lower Unit range between 24 and 16 Ma, which is consistent with a protracted presence of melt in the Lower Unit for about 8 Myr. The spread of the zircon U-Pb ages along the concordia from 24 to 16 Ma is interpreted as a succession of dissolution-crystallization cycles of 1–2 Myr, which can then be considered to represent the characteristic time of zircon revolution in a convection cell (Vanderhaeghe *et al.*, 2018). Granitic dikes intruding the Middle Unit form a network structurally connected to the migmatites of the Lower Unit (Vanderhaeghe, 2004). Dikes with a preserved magmatic fabric are oriented preferentially perpendicular or parallel to
the regional mineral and stretching lineation, which is consistent with intrusion in a context of pure flattening. These dikes cross-cut other dikes that are partially to totally transposed in the regional foliation, which suggests coeval intrusion and deformation. In the XZ plane of deformation (perpendicular to foliation and parallel to lineation), asymmetric boudinage of the dikes is consistent with a top to the NNE sense of shear and thus intrusion during regional extension. In the YZ plane of deformation, perpendicular to the NNE-SSW mineral and stretching lineation, asymmetric boudinage shows an opposite sense of shear on each side of the dome, top to the west on the western limb and top to the east on the eastern one. This transposition of granitic dikes rooted in the migmatites has been attributed to the formation of the first order dome during regional crustal extension (Vanderhaeghe, 2004). Zircon and monazite U-Pb ages range from 16 Ma, on an almost totally transposed dike, to 14 Ma on a less transposed one, which is consistent with the development of the Naxos first order dome over at least 2 Myr (Vanderhaeghe, 2004; Vanderhaeghe et al., 2018).

The petrological, structural and geochronological data presented above indicate a contrasted record between the Middle Unit and the Lower Unit in Naxos. Alltogether this record indicates that the Naxos domes were formed during the transition between crustal thickening and crustal extension. In particular, the dominant foliation of the Middle Unit, associated with F₁/F₂ fold interferences, is tilted and transposed into the syn-migmatitic foliation that delineates the nested dome structure of the migmatites of the Lower Unit. This structural record suggests that (i) dome formation postdates with F₁/F₂ fold interferences and (ii) the partially molten Lower Unit was mechanically decoupled from the Middle Unit during formation of the nested migmatite domes. These data serve to constrain the initial and boundary conditions of the numerical model presented in the next section on the one hand, and then to evaluate the characteristic length and time scales of Naxos's nested domes via the modeled scenario.

## 3 Numerical method

The numerical experiments are carried out with the open-source finite-volume code OpenFOAM. OpenFOAM is a C++ toolbox aimed at solving fluid mechanics problems. It includes a VOF method, which is an Eulerian fixed grid approach that tracks material-phases (fluids of distinct mechanical properties) interfaces, hence enabling it to compute large viscous deformations. The Navier-Stokes equations are solved together with the transport equations for the volume fraction of the fluid phases. As shown in Louis-Napoléon et al. (2020) this approach gives accurate results for standard Rayleigh-Taylor and Rayleigh-Bénard problems, and in a reasonable computational time (domain decomposition and the Message Passing Interface (MPI) libraries are used to increase the computational speed). We have developed our own solver (cf. Louis-Napoléon et al., 2020, 2022, 2024), which can be downloaded at: https://gitlab.com/AurelieLN/MultiMeltInterFoam; equations are recalled below. Details of the computational technique can be found in Louis-Napoléon et al. (2020, 2022, 2024) and will not be repeated here, for conciseness.

The temperature of the modeled crust is solved via the heat equation (Carslaw and Jaeger, 1959; England and Thompson, 1984):

$$ \rho_{ref} C_p \left[ \frac{\partial T}{\partial t} + U \cdot \nabla T \right] = - \nabla \cdot (\rho_{ref} C_p(T) \kappa(T) \nabla T) + H_r \quad (1) $$

where T is temperature (K), ρ is the local density (kg m⁻³), U the local velocity (m.s⁻¹), κ is the thermal diffusivity (m² s⁻¹), H is the radiogenic heat production (W m⁻³), and Cₚ is the heat capacity (J kg⁻¹ K⁻¹). Here we have not considered the influence of latent heat upon melting, but its influence has been discussed in our previous work (Louis-Napoléon et al., part II, 2024); we showed there that it can affect heat transfer by about 20%, hence slightly shift the flow regimes ranges. Yet it does not modify, in the present test case, the resulting model dynamics. The heat equation is coupled to the momentum equation to solve for the local velocity field U, accounting for pressure P and gravity g:

$$ \rho_{ref} \frac{\partial U}{\partial t} + \rho_{ref} U \cdot \nabla U = - \nabla P + \rho g + \nabla \cdot \left[ \mu (\nabla U + (\nabla U)^T) \right] \quad (2) $$

The viscosity of the ductile crust follows a creep power-law function of temperature (T) and shear strain rate ˙ε (Chen and Morgan, 1990):

$$ \mu = K_{eff}(T) \dot{e}^{\frac{1}{n}-1} \quad \text{with} \quad K_{eff} = \sum_{i=1}^{3} C_i \left[ 0.25 \times (0.75A_i)^{\frac{1}{n}} \times \exp\left(\frac{Q}{nRT}\right) \right] \quad (3) $$

where Kₐ₇₅ is the consistency (in kg m⁻¹ s⁻²+1/n) and corresponds to a dynamic viscosity when n=1. ˙ε has a minimal value set to 10⁻¹⁶ s⁻¹, R=8.314 J mol⁻¹ K⁻¹ is the universal gas constant, and Aᵢ is a material phase dependent pre-factor, with Cᵢ the volume fraction of each phase i. The numerical term [0.25 × (0.75)⁻¹] stems from adapting the uniaxial strain rate flow law deduced from lab experiments to an isotropic tensor stress-strain rate relation for incompressible material (Chen and Morgan, 1990, cf. their equations 1 to 5). Activation energy, Q = 1.54 × 10⁵ J mol⁻¹, and n = 2.3 are taken from mechanical experiments performed on quartz (Ranalli, 1995). The choice of this dominant composition stems from the very felsic composition displayed by Naxos's crustal material (see description of the geological context).

Here, we assume that partial melting results in contrasted density (Δρ up to about 300 kg/m³) and viscosity (Δμ ranging over 10⁵ Pa.s) of the white and black layers with respect to the ambient medium, which mimics melt percolation along grain boundaries and small veins followed by progressive gathering and accumulation into clusters or layers up to hectometer in size (e.g. Räss et al., 2019; Edmonds et al., 2019). This assumption presents the advantage to maintain a single flow law in the model formulation that, however, accounts for distinct density and viscosity parameters for each material phase. Note that in Louis-Napoléon et al. Part I (2022) and Part II (2024) we explicitly incorporated a melt dependent viscosity
law that is activated above the 30% melt threshold; the simpler assumption that we make here aims at showing that considering a simple and drastic viscosity drop upon crustal melting is actually sufficient to model the process of convection, diapirism and nested dome formation. We consider that this assumption is a good balance between i) a simple model with few parameters and ii) a model that still grasps the first order mechanical impact of partially melting rocks (Ganne et al., 2014; Vanderhaeghe et al., 2003). We adopt Boussinesq's approximation by accounting for variable density only in the gravity term of the momentum equation (2); density is indeed supposed to obey the temperature-dependent state equation:

$$ \rho = \rho_{ref} \times [1 - \alpha(T - T_{ref})] \\ \text{with } \rho_{ref} = \sum_{i=1}^{3} C_i \rho_{refi} \quad (4) $$

where $\alpha$ is the thermal expansion ($\alpha = 3 \times 10^{-5} \text{ K}^{-1}$), $T_{ref}$ is the temperature imposed at the top of the domain, and $\rho_{ref}$ is a local phase-dependent density which depends on $C_i$ (phase i volume fraction), and $\rho_{refi}$ is a reference density for each material phase (see definition below).

## 4 Numerical setup

The modeled crustal domain is 45 km thick and 50 km wide. It is subdivided into a 10 km thick rigid upper crust overlying a viscous crust made of layers with contrasting physical properties relative to the ambient medium, namely some layers are more dense and more viscous, whereas others are less dense and viscous (black versus white layers, respectively). These layers aim at representing the impact of hectometer-scale heterogeneities on the behavior of a partially molten crust. Instead of implementing the complex physics that drive sub-scale melt transfer through a crust that starts to melt by percolation through the effective porous media (e.g. Räss et al., 2019; Schmeling et al., 2019), we assume that upon crustal melting, this sub-scale melt transport forms a heterogeneous medium of characteristic size of several hundred meters (e.g. Duretz et al., 2019; Edmonds et al., 2019). In previous investigations, different sizes, shapes and concentration of inclusions have been tested (Martin and Nokes, 1989; Harada et al., 2012; Louis-Napoléon et al., 2024), showing that various flow regimes could develop, ranging from local segregation of the inclusions to their suspension within large scale convective flow. In Louis-Napoléon et al. (2022), we also showed that a specific regime allowed for simultaneous convection and layering of the buoyant inclusions, for appropriate characteristic Rayleigh numbers, inclusions sizes greater than 300 m and concentrations greater than 0.35. In that regime, spherical inclusions deform, disperse during the development of crustal-scale gravitational instabilities, but the buoyant ones can also re-aggregate to form domes and layers stacking at the top of the convective cells. In the present paper, we report the result of a model with initial layers rather than with initial spherical inclusions, which represents a more realistic geological setting for Naxos (Fig. 2). These layers represent more felsic (white) and more mafic (black) lithologies within a melting crust of intermediate composition, designated as the ambient medium.

<figure><img src="image_2.png" /><figcaption>Fig. 2. Initial set-up of the numerical model. The upper crust above 10 km depth is rigid and is not represented. From -10 km to -45 km depth, the behavior of the modeled crust follows a power-law temperature and strain rate dependent rheology. Heterogeneities are represented by layers that are more dense and viscous (black) and less dense and viscous (white) than the ambient medium.</figcaption></figure>

The ambient medium, the white layers and the black layers are denoted hereafter as material “phases” which are identified by volume fractions $C_1$, $C_2$ and $C_3$, respectively.

The initial temperature linearly increases from $T_c = 300 \text{ }^{\circ}\text{C}$ at -10 km depth to $T_h = 600 \text{ }^{\circ}\text{C}$ at -45 km depth. At the onset of the numerical experiment, the basal temperature is switched to $T_{h+} = 1000 \text{ }^{\circ}\text{C}$, simulating the thermal impact of slab break-off or delamination (Ueda et al., 2012; Vanderhaeghe and Duchêne, 2010).

The ambient medium and the layers are assigned different values of heat production, density and parameter A, which controls the viscosity (Tab. 1). White layers produce more heat and are less dense and less viscous than the ambient medium, mimicking felsic rocks. In contrast, black layers produce less heat and are denser and more viscous than the ambient medium, mimicking mafic rocks. A detailed parametric analysis of the model's sensitivity to the rheological formulation, the temperature evolution and the size and distribution of inclusions is provided in Louis-Napoléon et al. (2022, part I) and Louis-Napoléon et al. (part II, 2024).

## 5 Results of the numerical model

The time evolution of our reference numerical experiment is presented in Figures 3 and 4. Heat diffuses from the base of the domain across the modeled domain. At ca. 2 Myr, domains where the viscosity of the ambient medium decreases below ~ $10^{19} \text{ Pa.s}$
<table><caption>**Table 1.** Parameters used for the computation of the density, the parameter A entering into the computation of viscosity and heat production of the ambient medium, white (felsic lithologies) and black (mafic lithologies) layers. Note that the viscosity pre-factor A (cf. Eq. 3) decreases with increasing viscosity. For absolute values of the corresponding viscosities, please refer to subsequent model figures. Parameters taken from Turcotte and Schubert (2014).</caption><thead><tr><th></th><th>ρ<sub>ref</sub> (kg/m<sup>3</sup>)</th><th>A (Pa<sup>-2.3</sup>s<sup>-1</sup>)</th><th>H<sub>r</sub> (W/m<sup>3</sup>)</th></tr></thead><tbody><tr><td>Ambient medium</td><td>2750</td><td>10<sup>18</sup></td><td>1.75</td></tr><tr><td>White layers</td><td>2600</td><td>10<sup>20</sup></td><td>2</td></tr><tr><td>Black layers</td><td>2900</td><td>5.10<sup>15</sup></td><td>0.9</td></tr></tbody></table>

(corresponding to a temperature above ~ 640 °C), the gravitationally unstable layers start to move relative to the ambient medium. At the bottom of the modeled crust, white, less dense and less viscous material migrates upward whereas black, denser and more viscous material, moves downward. The accumulation of low-density pockets below the lowermost black layer is followed by the development of a composition-driven diapir at ca. 2.5 Myr. At the top of the diapir, the white layer is aggregated into pockets and the black layer is disaggregated. Similarly, the accumulation of low-density white pockets below the shallowest dense and viscous black layer triggers the development of a diapir at ca. 5 Myr. Thermally-driven convection starts at ca. 10 Ma, when the 640 °C isotherm has reached – 25 km depth and about half of the modeled crust has a viscosity below ~ 10¹⁹ Pa·s. At this stage, the convection cell is about 15 km in diameter. As temperature continues to rise in the crust, the diameter of the convection cells, marked by deflection of the isotherms, increases to about 25 km at ca. 20 Myr. As shown in Fig. 4, markers of material, initially located at distinct depths, are entrained in the convection cells with revolving cycles of 2.5 Myr on average. The velocity of these markers ranges from 2 cm/yr to 4 cm/yr and their temperature oscillates around 800 °C +/- 100 °C. The less dense material clustered above the convection cells forms dome structures of diameter ~ 5 km. The black and denser material, in turn, has settled down at the base of the convective cell, and forms a ~ 5 km thick homogeneous layer.

# 6 Discussion

In this section, we discuss the pertinence of the assumptions and boundary conditions of the model and their implications in the context of the thermo-mechanical evolution of the Aegean domain. First, let us note that the characteristic length and time scales of the modeled nested domes are comparable with the geological record of Naxos's nested domes.

The basal heating condition, marked by an instantaneous switch from 600 °C to 1000 °C, simulates heat advection from the asthenosphere that would replace the lithospheric mantle. This situation applies to the Aegean domain marked by convergence and slab retreat during the Cenozoic (Jolivet and Brun, 2010; Ring et al., 2010; Spakman et al., 1988; Vanderhaeghe et al., 2007). The persistence of this high temperature over 20 Myr is justified by the thermal relaxation time of the lithosphere, which may last several tens of Myr in such a context of slab retreat or delamination (Ueda et al., 2012). In other words, in such a

situation, heat advection dominates over diffusion for several tens of Myr, which justifies the assumption of constant basal temperature in our model. Further discussion on this thermal boundary condition can be found in Louis-Napoléon *et al.* (2024).

No-slip boundary conditions are applied at the upper (– 10 km depth) and lower boundaries (-45 km depth) of the model domain, respectively, while the lateral borders are assigned periodic boundary conditions (as in Louis Napoléon *et al.*, 2022). This mimics a transitional stage during orogenic evolution, which might correspond to the development of an orogenic plateau, after tectonic accretion and before gravitational collapse (Vanderhaeghe, 2009, 2012). The periodic lateral boundary conditions assumption implies that the lateral extent of basal heating and geometric distribution of the heterogeneous layers extends over a much broader area than just the model domain. The extent of partial melting in the Aegean domain might be apprehended by the distribution of HT/LP metamorphic rocks and Cenozoic granitic plutons that spread over the central part of the Cyclades (Vanderhaeghe and Teyssier, 2001; Rabillard *et al.*, 2018; Vanderhaeghe *et al.*, 2007; Jolivet *et al.*, 2021).

Experimental petrology and thermodynamic modeling indicate that the onset of partial melting of crustal rocks in presence of water occurs at temperatures ranging from 600 to 650 °C (Collins *et al.*, 2021; Weinberg and Hasalová, 2015). Partial melting by destabilization of micas takes place at temperatures exceeding 650 °C (Gardien *et al.*, 1995; Patino Douce and Johnston, 1991; Vielzeuf and Holloway, 1988) while destabilization of amphiboles happens at about 700 °C (Palin *et al.*, 2016; Palin *et al.*, 2016; Rapp *et al.*, 1991). These investigations show that depending on composition, melt fractions reach 20-40 % at about 900 °C (Patino Douce and Johnston, 1991; Thompson and Connolly, 1995). Such melt fractions-temperature relationships are consistent with those identified in the field at Naxos (Kruckenberg *et al.*, 2011; Lamont *et al.*, 2023; Peillod *et al.*, 2017).

The rheology of partially molten rocks and magmas is known to be controlled by the relative proportions of melt and solids and is characterized by two distinct thresholds (Arzi, 1978; Rosenberg, 2001; Rutter *et al.*, 2006; van der Molen and Paterson, 1979; Vanderhaeghe, 2009; Vigneresse *et al.*, 1996). The onset of partial melting, with a melt fraction of only a few percent, is marked by a strength decrease of about two orders of magnitude, while the transition from partially molten rock to magma occurs at about 30 % melt and is marked by an even more drastic strength loss (Arzi, 1978; Rutter *et al.*, 2006; van der Molen and Paterson, 1979). On the other hand, magma viscosity increases drastically at crystal proportions of about 70 % (Roscoe, 1952). In the model presented here, the viscosity of the ambient medium ranges from 10²³ at ca. 650°C to 10¹⁶ Pa·s at 900°C, which covers the range of values expected in a partially molten crust such as in Naxos. While in Louis-Napoléon *et al.* Part I (2022) and part II (2024) we have explicitly incorporated a melt dependent viscosity law above the 30% melt threshold, here we have chosen on purpose a simpler parametrization of viscosity based on the standard power law stress-strain relationship determined from lab experiments on rocks. We reckon that it is useful to consider more realistic melt dependent rheologies, and several studies have actually implemented coupled
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** 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^{16}$ to $10^{23}$ 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^{19}$ 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</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Modeled displacement (top), thermal evolution (middle) and velocity (bottom) of three material points, located initially at depths equal to -44 km (red), -40 km (blue), and -38 km (green), during the evolution of the 3D numerical experiment. The first accelerated motion is recorded after ca. 2 Myr of heating and corresponds to downward motion of the blue marker, concomitant with upward motion of the green marker, which corresponds to the onset of the buoyancy-driven segregation. Just before 2.5 Myr, the onset of composition-driven diapirism is marked by rapid motion of all markers, some going up (red and blue points), and others going down (green) at a velocity up to 15 cm/yr. A second rapid, gravitational destabilization of the layers occurs at 5 Myr. From this time onward, the three markers record cycles of temperature increase and decrease from 650 to 900 °C with a period of 1 to 3 Myr, and at a velocity of 1 to 7 cm/yr.</figcaption></figure>

formulations for two phase flow media (Keller et al., 2000; Schmeling et al., 2019, 2023) for other modeling purposes. However, one must also acknowledge that such implementations require many additional parameters and are non-unique, with rock lab experiments showing a wide variety of behaviors depending on composition, P, T, and apparatus conditions (Ganzhorn et al., 2016; Rosenberg, 2001; Rosenberg and Handy, 2005; Zhou et al., 2017). Therefore, we advocate that the simple parametrization proposed here based on layers with contrasted physical properties evidences the key influence of a geometrically heterogeneous distribution of density and viscosity on the formation of nested migmatite domes.

The density of crustal metamorphic rocks generally increases with depth (Bousquet et al., 1997; Hacker et al., 2003; Tassara, 2006). However, partial melting is marked by a decrease in density of about 10 % (e.g. Tassara, 2006). The partially molten crust is thus gravitationally unstable, and this

instability is further amplified by melt/solid segregation (Vanderhaeghe, 2009).

The consideration of layers in our model allows us to investigate the influence of heterogeneities within the crust that are inherited from sedimentary and magmatic processes and further enhanced by deformation, metamorphism and partial melting. While the thickness of these layers together with the resolution of the model do not allow to capture melt segregation processes at smaller scales, their 600 m thickness corresponds to the characteristic size that has previously been identified for melt segregation by compaction in felsic crust, potentially leading to the formation of networks of concordant-discordant veins (Brown et al., 1995; Edmonds et al., 2019; Räss et al., 2019; Vanderhaeghe, 2001). Louis Napoléon et al. (2022) showed with a series of tests that these heterogeneities tend to deform, disaggregate and then re-aggregate during the development of gravitational instabilities, so that their initial size and shape do not have a major impact on the model's
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Schematic representation of the formation of Naxos’s nested domes inspired by the numerical experiments. Partial melting of the crust allows upward migration of buoyant material (melt, magma, felsic partially molten layers represented in white) and accumulation of dense material (residual solids, cumulates, mafic partially molten layers in dark brown). The upward migration of the melting front is increasing the thickness of the convective domain with time, increasing the Rayleigh number and permitting larger convection cells to form and causing an additional decrease in viscosity and density. The first order dome is controlled by the size of the convection cells (drawn in pink) and the second order domes form within the low-density material aggregated at the top of the convection cells (drawn in white).</figcaption></figure>

evolution; the latter essentially shifts the onset of local segregation and/or regional convection by up to two million years (10% of the modeled basal heating time).

The model presented in this paper was specially designed to test the pertinence of the proposition that Naxos’s nested domes result from the development of gravitational instabilities. To emphasize the effects of buoyancy, we chose not to apply any lateral boundary conditions that would reproduce the influence of plate tectonics, nor vertical boundary conditions that would reproduce either mantle dynamics or top surface erosion/sedimentation processes. As a consequence, domes formed during the numerical experiment display a vertical axis of symmetry. This contrasts with the elliptical shape of Naxos’s main dome, which is thus interpreted to reflect the contribution of plate tectonics marked by North-East extension in the Aegean domain. Most of the exhumation of the former partially molten crust occurred on Naxos after ca. 12 Ma, namely during crustal extension, and has been attributed to isostasy-driven low-viscosity flow in the space left opened by the extension (Rey et al., 2011). However, considering the structural characteristics of the migmatites described above, the nested nature of the second order domes and their axisymmetric shape attest to the key role of buoyancy-driven forces (Vanderhaeghe et al., 2018). Moreover, the U-Pb geochronological zircon record from the migmatites coring the Naxos dome, with ages as old as ca. 24 Ma, indicates that these gravitational instabilities formed at least 8 Myr before the onset of Miocene regional extension (Vanderhaeghe et al., 2018).

## 7 Conclusion

Despite some limitations discussed above, the numerical experiment presented in this paper confirms that it is possible

to reproduce the structural and geochronological record of the migmatites coring Naxos domes with gravitational instabilities (Fig. 5), combining segregation, convection and diapirism. Namely, the modeled dome, which has a diameter of about 20 km, would reflect crustal-scale convection at a viscosity lower than $10^{18}$ Pa.s (consistent with partially molten rocks). The modeled revolving cycles of about 2.5 Myr are similar to the cyclic geochronological record of sampled zircon grains in the migmatites. In the model, the second order domes of ca. 5 km in diameter result from diapiric instabilities that develop both during and after the segregation of the light material above the convective cells. These domes are nested in the ~20 km wide broad dome.

Thermomechanical modeling of the formation of Naxos migmatite-cored domes demonstrates that the development of gravitational instabilities within a low viscosity partially molten crust is an efficient mechanism to redistribute heterogeneities at the crustal scale, resulting in the formation of nested domes within a main dome of several tens of kilometers in size. Hence, these gravitational instabilities contribute to crustal-scale differentiation.

## Acknowledgments

This paper is part of Aurélie Louis-Napoléon’s PhD thesis. Aurélie Louis-Napoléon has been supported by a PhD fellowship from the French Ministry of Higher Education and Research and by research credits from the CNRS-INSU Syster program. This work was granted access to the HPC resources of CALMIP supercomputing center (https://www.calmip.univ-toulouse.fr/) under the allocation reference P1525-2023. We would like to thank Alexander Cruden, Andrew Zuza and Laetita Le Pourhiet their constructive reviews, which resulted in significant improvement of the paper.
References

Arzi AA. 1978. Critical phenomena in the rheology of partially melted rocks. *Tectonophysics* 44: 173-184.

Avigad D. 1998. High-pressure metamorphism and cooling on SE Naxos (Cyclades, Greece). *Eur J Mineral* 10: 1309-1320.

Avigad D, Ziv A, Garfunkel Z. 2001. Ductile and brittle shortening, extension-parallel folds and maintenance of crustal thickness in the central Aegean (Cyclades, Greece). *Tectonics* 20: 277-287.

Babeyko AY, Sobolev S, Trumbull R, Oncken O, Lavier L. 2002. Numerical models of crustal scale convection and partial melting beneath the Altiplano-Puna plateau. *Earth Planet Sci Lett* 199: 373-388.

Bousquet R, Goffé, B., Henry P, Le Pichon X, Chopin C. 1997. Kinematic, thermal and petrological model of the Central Alps: lepontine metamorphism in the upper crust and eclogitisation of the lower crust. *Tectonophysics* 273: 105-127.

Brichau S, Ring U, Ketcham RA, Carter A, Stockli D, Brunel M. 2006. Constraining the long-term evolution of the slip rate for a major extensional fault system in the central Aegean, Greece, using thermochronology. *Earth Planet Sci Lett* 241: 293-306.

Brown M. 2001. Orogeny, migmatites and leucogranites: a review. *J Earth Syst Sci* 110: 313-336.

Brown M, Averkin Y.A, McLellan E.L, Sawyer E.W. 1995. Melt segregation in migmatites. *J. Geophys. Res. Solid Earth* 100, 15655-15679. https://doi.org/10.1029/95JB00517

Brun J-P. 1983. L'origine des domes gneissiques; modèles et tests. *Bull. Société Géologique Fr.* S7-XXV: 219-228.

Brun J-P., Gapais D, Le Theoff B. 1981. The mantled gneiss domes of Kuopio (Finland): interfering diapirs. *Tectonophysics* 74: 283-304.

Brun J-P., Sokoutis D, Van Den Driessche F J. 1994. Analogue modeling of detachment fault systems and core complexes. *Geology* 22: 319-322.

Buck WR. 1991. Modes of continental lithospheric extension. *J Geophys Res Solid Earth* 96: 20161-20178.

Buick IS. 1991. The late Alpine evolution of an extensional shear zone, Naxos, Greece. *J Geol Soc* 148: 93-103.

Buick IS, Holland TJB. 1989. The P-T path associated with crustal extension, Naxos, Cyclades, Greece. *Geol Soc Lond Spec Publ* 43: 365-369.

Burg J-P., Kaus BJP, Podladchikov YY. 2004. Dome structures in collision orogens: Mechanical investigation of the gravity/ compression interplay, in: Gneiss Domes in Orogeny. *Geological Society of America*. https://doi.org/10.1130/0-8137-2380-9.47

Burg J-P, Podladchikov Y. 2000. From buckling to asymmetric folding of the continental lithosphere: numerical modelling and application to the Himalayan syntaxes. *Geol Soc Lond Spec Publ* 170: 219-236.

Carslaw HS, Jaeger JC. 1959. *Conduction of heat in solids*, 2nd ed. Oxford: Clarendon Press.

Chen Y, Morgan WJ. 1990. A nonlinear rheology model for mid-ocean ridge axis topography. *J Geophys Res* 95: 17583.

Collins WJ. 1989. Polydiapirism of the Archean Mount Edgar Batholith, Pilbara Block, Western Australia. *Precambrian Res* 43: 41-62.

Collins WJ, Murphy JB, Blereau E, Huang H-Q. 2021. Water availability controls crustal melting temperatures. *Lithos* 402-403: 106351.

Coney PJ, Harms TA. 1984. Cordilleran metamorphic core complexes: Cenozoic extensional relics of Mesozoic compression. *Geology* 12: 550.

Cruden AR, Koyi H, Schmeling H. 1995. Diapiric basal entrainment of mafic into felsic magma. *Earth Planet Sci Lett* 131: 321-340.

Davis GH. 1983. Shear-zone model for the origin of metamorphic core complexes. *Geology* 11: 342.

Dercourt J, Zonenshain L.P., Ricou L.-E., Kazmin V.G, Le Pichon X, Knipper A.L, Grandjacquet C, Sbortshikov I.M, Geyssant J, Lepvrier C, Pechersky D.H, Boulin J, Sibuet J.-C, Savostin L.A, Sorokhtin O, Westphal M, Bazhenov M.L, Lauer J.P, Biju-Duval B. 1986. Geological evolution of the tethys belt from the atlantic to the pamirs since the LIAS. *Tectonophysics* 123, 241-315. https://doi.org/10.1016/0040-1951(86)90199-X

Dewey J.F., Şengör A.M.C. 1979. Aegean and surrounding regions: Complex multiplate and continuum tectonics in a convergent zone. *Geol. Soc. Am. Bull.* 90, 84. https://doi.org/10.1130/0016-7606 (1979)90<84:AASRCM>2.0.CO;2

Duchêne S, Aïssa R, Vanderhaeghe O. 2006. Pressure-Temperature-time Evolution of Metamorphic Rocks from Naxos (Cyclades, Greece): constraints from Thermobarometry and Rb/Sr dating. *Geodin Acta* 19: 301-321.

Duretz T, Räss L, Podladchikov Y, Schmalholz S. 2019. Resolving thermomechanical coupling in two and three dimensions: spontaneous strain localization owing to shear heating. *Geophys J Int* 216: 365-379.

Edmonds M, Cashman KV, Holness M, Jackson M. 2019. Architecture and dynamics of magma reservoirs. *Philos Trans R Soc Math Phys Eng Sci* 377. 20180298.

England PC, Thompson AB. 1984. Pressure—temperature—time paths of regional metamorphism I. Heat transfer during the evolution of regions of thickened continental crust. *J Petrol* 25: 894-928.

Eskola PE. 1948. The problem of mantled gneiss domes. *Q J Geol Soc* 104: 461-476.

Fytikas M, Innocenti F, Manetti P, Peccerillo A, Mazzuoli R, Villari L. 1984. Tertiary to quaternary evolution of volcanism in the Aegean region. *Geol Soc Lond Spec Publ* 17: 687-699.

Ganne J, Gerbault M, Block S. 2014. Thermo-mechanical modeling of lower crust exhumation—Constraints from the metamorphic record of the Palaeoproterozoic Eburnean orogeny, West African Craton. *Precambrian Res* 243: 88-109.

Ganzhorn AC, Trap P, Arbaret L, Champallier R, Fauconnier J, Labrousse L, Prouteau G. 2016. Impact of gneissic layering and localized incipient melting upon melt flow during experimental deformation of migmatites. *J Struct Geol* 85: 68-84.

Gardien V, Thompson AB, Grujic D, Ulmer P. 1995. Experimental melting of biotite+ plagioclase+ quartz±muscovite assemblages and implications for crustal melting. *J Geophys Res Solid Earth* 100: 15581-15591.

Gautier P, Brun J-P., Jolivet L. 1993. Structure and kinematics of Upper Cenozoic extensional detachment on Naxos and Paros (Cyclades Islands, Greece). *Tectonics* 12: 1180-1194.

Gautier P, Brun, J-P., Moriceau R, Sokoutis D, Martinod J, Jolivet L. 1999. Timing, kinematics and cause of Aegean extension: a scenario based on a comparison with simple analogue experiments. *Tectonophysics* 315: 31-72.

Hacker BR, Abers GA, Peacock SM. 2003. Subduction factory 1. Theoretical mineralogy, densities, seismic wave speeds, and H₂O contents: subduction zone mineralogy and physical properties. *J Geophys Res Solid Earth* 108. https://doi.org/10.1029/2001JB001127
Harada S, Mitsui T, Sato K. 2012. Particle-like and fluid-like settling of a stratified suspension. *Eur Phys J E* 35. https://doi.org/10.1140/epje/i2012-12001-6

Jansen JB. 1973. *Geological Map of Naxos*.

Jansen JBH, Schuiling RD. 1976. Metamorphism on Naxos; petrology and geothermal gradients. *Am J Sci* 276: 1225-1253.

Jolivet L, Brun J-P. 2010. Cenozoic geodynamic evolution of the Aegean. *Int J Earth Sci* 99: 109-138.

Jolivet L, Arbaret L, Le Pourhiet L, Cheval-Garabédian F, Roche V, Rabillard A, Labrousse L. 2021. Interactions of plutons and detachments: a comparison of Aegean and *Tyrrhenian granitoids*. *Solid Earth* 12: 1357-1388.

Keay S, Lister G, Buick I. 2001. The timing of partial melting, Barrovian metamorphism and granite intrusion in the Naxos metamorphic core complex, Cyclades, Aegean Sea, Greece. *Tectonophysics* 342: 275-312.

Keller AA, Blunt MJ, Roberts PV. 2000. Behavior of nonaqueous phase liquids in fractured porous media under two-phase flow conditions. *Transp Porous Media* 38: 189-203.

Kruckenberg SC, Ferré EC, Teyssier C, Vanderhaeghe O, Whitney DL, Seaton NCA, Skord JA. 2010. Viscoplastic flow in migmatites deduced from fabric anisotropy: an example from the Naxos dome, Greece. *J Geophys Res* 115. https://doi.org/10.1029/2009JB007012

Kruckenberg SC, Vanderhaeghe O, Ferré EC, Teyssier C, Whitney DL. 2011. Flow of partially molten crust and the internal dynamics of a migmatite dome, Naxos, Greece: internal dynamics of the Naxos dome. *Tectonics* 30: n/a-n/a.

Lamont TN, Searle MP, Waters DJ, Roberts NMW, Palin RM, Smye A, Dyck B, Gopon P, Weller OM, St-Onge MR. 2020. Compressional origin of the Naxos metamorphic core complex, Greece: structure, petrography, and thermobarometry. *GSA Bull* 132: 149-197.

Lamont TN, Smye AJ, Roberts, N.M.W., Searle MP, Waters DJ, White RW. 2023. Constraints on the thermal evolution of metamorphic core complexes from the timing of high-pressure metamorphism on Naxos, Greece. *GSA Bull*. https://doi.org/10.1130/B36332.1

Le Pourhiet L, Huet B, May DA, Labrousse L, Jolivet L. 2012. Kinematic interpretation of the 3D shapes of metamorphic core complexes. *Geochem Geophys Geosyst* 13. https://doi.org/10.1029/2012GC004271.

Linnros H, Hansman R, Ring U. 2019. The 3D geometry of the Naxos detachment fault and the three-dimensional tectonic architecture of the Naxos metamorphic core complex, Aegean Sea, Greece. *Int J Earth Sci* 108: 287-300.

Lister GS, Banga G, Feenstra A. 1984. Metamorphic core complexes of Cordilleran type in the Cyclades, Aegean Sea, Greece. *Geology* 12: 221.

Louis-Napoléon A, Gerbault M, Bonometti T, Thieulot C, Martin R, Vanderhaeghe O. 2020. 3-D numerical modelling of crustal polydiapirs with volume-of-fluid methods. *Geophys J Int* 222: 474-506.

Louis-Napoléon A, Bonometti T, Gerbault M, Martin R, Vanderhaeghe O. 2022. Models of convection and segregation in heterogeneous partially molten crustal roots with a VOF method – I: flow regimes, *Geophys J Int* 229: 2047-2080.

Louis-Napoléon A, Gerbault M, Bonometti T, Vanderhaeghe O, Martin R, Maury N. 2024. Convection and segregation in heterogeneous orogenic crust with a VOF method – II: how to form migmatite domes. *Geophys J Int* 236: 207-232.

Martin D, Nokes R. 1989. A fluid-dynamical study of crystal settling in convecting magmas. *J Petrol* 30: 1471-1500.

Martin L, Duchêne S, Deloule E, Vanderhaeghe O. 2006. The isotopic composition of zircon and garnet: a record of the metamorphic history of Naxos, Greece. *Lithos* 87: 174-192.

Martin LAJ, Duchêne S, Deloule E, Vanderhaeghe O. 2008. Mobility of trace elements and oxygen in zircon during metamorphism: consequences for geochemical tracing. *Earth Planet Sci Lett* 267: 161-174.

Myers JS, Watkins KP. 1985. Origin of granite-greenstone patterns, Yilgarn Block, Western Australia. *Geology* 13: 778.

Palin RM, White RW, Green ECR. 2016. Partial melting of metabasic rocks and the generation of tonalitic-trondhjemitic-granodioritic (TTG) crust in the Archaean: constraints from phase equilibrium modelling. *Precambrian Res* 287: 73-90.

Palin RM, White RW, Green ECR, Diener JFA, Powell R, Holland TJB. 2016. High-grade metamorphism and partial melting of basic and intermediate rocks. *J Metamorph Geol* 34: 871-892.

Patino Douce AE, Johnston AD. 1991. Phase equilibria and melt productivity in the pelitic system: implications for the origin of peraluminous granitoids and aluminous granulites. *Contrib Mineral Petrol* 107: 202-218.

Peillod A, Majka J, Ring U, Drüppel K, Patten C, Karlsson A, Włodek A, Fehler E. 2021a. Differences in decompression of a high-pressure unit: A case study from the Cycladic Blueschist Unit on Naxos Island, Greece. *Lithos* 386-387: 106043.

Peillod A, Ring U, Glodny J, Skelton A. 2017. An Eocene/Oligocene blueschist-/greenschist facies *P-T* loop from the Cycladic Blueschist Unit on Naxos Island, Greece: deformation-related re-equilibration vs. thermal relaxation. *J Metamorph Geol* 35: 805-830.

Peillod A, Fehler E, Ring U. 2021b. Quo vadis Zeus: is there a Zas shear zone on Naxos Island, Aegean Sea, Greece? A review of metamorphic history and new kinematic data. *J Geol Soc* 178. https://doi.org/10.1144/jgs2020-217

Petford N, Cruden AR, McCaffrey KJW, Vigneresse J-L. 2000. Granite magma formation, transport and emplacement in the Earth's crust. *Nature* 408: 669-673.

Porada H, Berhorst V. 2000. Towards a new understanding of the Neoproterozoic-early palæozoic Lufilian and northern Zambezi belts in Zambia and the Democratic Republic of Congo. *J Afr Earth Sci* 30: 727-771.

Rabillard A, Jolivet L, Arbaret L, Bessière E, Laurent V, Menant A, Augier R, Beaudoin A. 2018. Syntensional granitoids and detachment systems within cycladic metamorphic core complexes (Aegean Sea, Greece): toward a regional tectonomagmatic model. *Tectonics* 37. https://doi.org/10.1029/2017TC004697

Ramberg H. 1981. The role of gravity in orogenic belts. *Geol Soc Lond Spec Publ* 9: 125-140.

Ranalli G. 1995. *Rheology of the Earth*, Springer, ed.

Rapp RP, Watson EB, Miller CF. 1991. Partial melting of amphibolite/eclogite and the origin of Archean trondhjemites and tonalites. *Precambrian Res.* 51: 1-25.

Räss L, Duretz T, Podladchikov YY. 2019. Resolving hydromechanical coupling in two and three dimensions: spontaneous channelling of porous fluids owing to decompaction weakening. *Geophys J Int* 218: 1591-1616.

Rey PF, Teyssier C, Kruckenberg SC, Whitney DL. 2011. Viscous collision in channel explains double domes in metamorphic core complexes. *Geology* 39: 387-390.
Riel N, Mercier J, Weinberg R. 2016. Convection in a partially molten metasedimentary crust? Insights from the El Oro complex (Ecuador). *Geology* 44: 31–34.

Ring U, Glodny J, Will T, Thomson S. 2010. The hellenic subduction system: high-pressure metamorphism, exhumation, normal faulting, and large-scale extension. *Annu Rev Earth Planet Sci* 38: 45–76.

Roscoe R. 1952. The viscosity of suspensions of rigid spheres. *Br J Appl Phys* 3: 267–269.

Rosenberg CL. 2001. Deformation of partially molten granite: a review and comparison of experimental and natural case studies. *Int J Earth Sci* 90: 60–76.

Rosenberg CL, Handy MR. 2005. Experimental deformation of partially melted granite revisited: implications for the continental crust. *J Metamorph Geol* 23: 19–28.

Rutter EH, Brodie KH, Irving DH. 2006. Flow of synthetic, wet, partially molten "granite" under undrained conditions: an experimental study: FLOW OF PARTIALLY MOLTEN "GRANITE." *J Geophys Res Solid Earth* 111: n/a-n/a. https://doi.org/10.1029/2005JB004257

Sawyer EW. 1994. Melt segregation in the continental crust. *Geology* 22: 1019–1022.

Schenker FL, Gerya T, Burg J-P. 2012. Bimodal behavior of extended continental lithosphere: Modeling insight and application to thermal history of migmatitic core complexes. *Tectonophysics* 579: 88–103.

Schmeling H, Marquart G, Weinberg R, Kumaravel P. 2023. Dynamic two-phase flow modeling of melt segregation in continental crust: batholith emplacement versus crustal convection, with implications for magmatism in thickened plateaus. *Geochem Geophys Geosyst* 24. https://doi.org/10.1029/2023GC010860

Schmeling H, Marquart G, Weinberg R, Wallner H. 2019. Modelling melting and melt segregation by two-phase flow: new insights into the dynamics of magmatic systems in the continental crust. *Geophys J Int* 217: 422–450.

Seward D, Vanderhaeghe O, Siebenaller L, Thomson S, Hibsch C, Zingg A, Holzner P, Ring U, Duchêne S. 2009. Cenozoic tectonic evolution of Naxos Island through a multi-faceted approach of fission-track analysis. *Geol Soc Lond Spec Publ* 321: 179–196.

Spakman W, Wortel, M.J.R., Vlaar NJ. 1988. The Hellenic Subduction Zone: a tomographic image and its geodynamic implications. *Geophys Res Lett* 15: 60–63.

Talbot CJ. 1979. Infrastructural migmatitic upwelling in East Greenland interpreted as thermal convective structures. *Precambrian Res* 8: 77–93.

Tassara A. 2006. Factors controlling the crustal density structure underneath active continental margins with implications for their evolution: CONTINENTAL MARGIN CRUSTAL DENSITY. *Geochem Geophys Geosyst* 7: n/a-n/a.

Thompson AB, Connolly JA. 1995. Melting of the continental crust: some thermal and petrological constraints on anatexis in continental collision zones and other tectonic settings. *J Geophys Res Solid Earth* 100: 15565–15579.

Tirel C, Gueydan F, Tiberi C, Brun J-P. 2004. Aegean crustal thickness inferred from gravity inversion. Geodynamical implications. *Earth Planet Sci Lett* 228: 267–280.

Turcotte D, Schubert G. 2014. *Geodynamics*, 3rd ed. Cambridge University Press.

Ueda K, Gerya TV, Burg J-P. 2012. Delamination in collisional orogens: Thermomechanical modeling: DELAMINATION IN COLLISIONAL OROGENS. *J Geophys Res Solid Earth* 117: n/a–n/a.

Urai JL, Schuiling RD, Jansen JB. 1990. Alpine deformation on Naxos (Greece), in: *Deformation Mechanisms, Rheology and Tectonics*, Geological Society Special Publication. Knipe, R. J. & Rutter, E. H., pp. 509–522.

van der Molen F I, Paterson MS. 1979. Experimental deformation of partially-melted granite. *Contrib Mineral Petrol* 70: 299–318.

Van Kranendonk MJ, Collins WJ, Hickman A, Pawley MJ. 2004. Critical tests of vertical vs. horizontal tectonic models for the Archaean East Pilbara Granite-Greenstone Terrane, Pilbara Craton, Western Australia. *Precambrian Res* 131: 173–211.

Vanderhaeghe O. 2001. Melt segregation, pervasive melt migration and magma mobility in the continental crust: the structural record from pores to orogens. Phys. Chem. Earth Part Solid Earth Geod. 26, 213–223. https://doi.org/10.1016/S1464-1895 (01)00048-5

Vanderhaeghe O. 2004. Structural development of the Naxos migmatite dome, in: Gneiss Domes in Orogeny. *Geological Society of America*. https://doi.org/10.1130/0-8137-2380-9.211

Vanderhaeghe O. 2009. Migmatites, granites and orogeny: Flow modes of partially-molten rocks and magmas associated with melt/solid segregation in orogenic belts. *Tectonophysics* 477: 119–134.

Vanderhaeghe O, Duchêne S. 2010. Crustal-scale mass transfer, geotherm and topography at convergent plate boundaries: crustal dynamics at convergent plate boundaries. *Terra Nova* 22: 315–323.

Vanderhaeghe O. 2012. The thermal–mechanical evolution of crustal orogenic belts at convergent plate boundaries: A reappraisal of the orogenic cycle. J. Geodyn. 56–57, 124–145. https://doi.org/10.1016/j.jog.2011.10.004

Vanderhaeghe O, Hibsch C, Siebenaller L, Duchêne S, de St Blanquat M, Kruckenberg S, Fotiadis A, Martin L. 2007. Penrose Conference – Extending a Continent – Naxos Field Guide. *J. Virtual Explor*. 27. https://doi.org/10.3809/jvirtex.2007.00175

Vanderhaeghe O, Kruckenberg SC, Gerbault M, Martin L, Duchêne S, Deloule E. 2018. Crustal-scale convection and diapiric upwelling of a partially molten orogenic root (Naxos dome, Greece). *Tectonophysics* 746: 459–469.

Vanderhaeghe O, Medvedev S, Fullsack P, Beaumont C, Jamieson RA. 2003. Evolution of orogenic wedges and continental plateaux: insights from crustal thermal-mechanical models overlying subducting mantle lithosphere. *Geophys J Int* 153: 27–51.

Vanderhaeghe O, Teyssier C. 2001. Partial melting and flow of orogens. *Tectonophysics* 342: 451–472.

Vielzeuf D, Holloway JR. 1988. Experimental determination of the fluid-absent melting relations in the pelitic system. *Contrib Mineral Petrol* 98: 257–276.

Vigneresse JL, Barbey P, Cuney M. 1996. Rheological transitions during partial melting and crystallization with application to felsic magma segregation and transfer. *J Petrol* 37: 1579–1600.

Weinberg RF. 1997. Diapir-driven crustal convection: decompression melting, renewal of the magma source and the origin of nested plutons. *Tectonophysics* 271: 217–229.

Weinberg RF, Hasalová, P. 2015. Water-fluxed melting of the continental crust: a review. *Lithos* 212-215: 158–188.

Weinberg RF, Schmeling H. 1992. Polydiapirs: multiwavelength gravity structures. *J Struct Geol* 14: 425–436.

Whitney DL, Teyssier C, Vanderhaeghe O. 2004. Gneiss domes and crustal flow. *Gneiss Domes Orogeny* 380: 15.
Wijbrans JR, McDougall I. 1988. Metamorphic evolution of the Attic Cycladic Metamorphic Belt on Naxos (Cyclades, Greece) utilizing 40Ar/39Ar age spectrum measurements. *J Metamorph Geol* 6: 571–594.

Yin A. 2004. Gneiss domes and gneiss dome systems, in: Gneiss Domes in Orogeny. *Geological Society of America*. https://doi.org/10.1130/0-8137-2380-9.1

Zhou Y, Zhang H, Yao W, Dang J, He C. 2017. An experimental study on creep of partially molten granulite under high temperature and wet conditions. *J Asian Earth Sci* 139: 15–29.

Zuza A, Cao W. 2023. Metamorphic core complex dichotomy in the North American Cordillera explained by Buoyant upwelling in variably thick crust. *GSA Today* 33: 4–11.

**Cite this article as:** Louis-Napoléon A, Vanderhaeghe O, Gerbault M, Martin R, Bonometti T. 2024. Formation of the Naxos nested domes and crustal differentiation by convection and diapirism, *BSGF - Earth Sciences Bulletin* **195**: 21.

Appendix

Here we recall the main results of our first analytical estimate of the conditions for crustal convection from Vanderhaeghe et al. (2018), Fig. A1, and of the regimes diagram obtained as a function of 2 characteristic Rayleigh numbers for partially melting crust, from Louis-Napoléon et al. (2022), Fig. A2.

Louis-Napoléon et al. (2022) identified, with numerical models, four different regimes of gravitational instabilities in a crust heated from below and containing heterogeneous inclusions (some more dense and viscous, other less dense and viscous than an ambient medium) that are evaluated on the basis of two characteristic Rayleigh numbers:

$$
\begin{align}
R_{aUM} &= \frac{2\left(\frac{H_T}{2}\right)^2}{\kappa_{UM}} \left( \frac{\rho\alpha g \Delta T_{UM} \left(\frac{H_T}{2}\right)}{2\kappa_{eff}^{UM}} \right)^n, \\
R_{aPM} &= \frac{\rho\alpha g \Delta T_{PM} \left(\frac{H_T}{2}\right)^3}{\kappa_{PM} \kappa_{eff}^{PM}} \tag{5}
\end{align}
$$

with $H_T$ half of the crust's height, $\rho$ density, $\alpha$ the thermal expansion coefficient, $\Delta T_{UM}$ and $\Delta T_{PM}$ the temperature gradient across the unmolten (UM) and partially molten (PM) crustal domains, $\kappa_{UM}$ and $\kappa_{PM}$ the corresponding thermal conductivities, and $K_{eff}^{UP}$ and $\kappa_{eff}^{PM}$ the corresponding consistency (equivalent to the inverse of a viscosity depending on its power law exponent).

<figure><img src="image_6.png" /><figcaption>**Fig. A1.** Analytical estimate of the conditions on crustal thickness (H) and corresponding average viscosity for crustal convection from Vanderhaeghe et al., (2018), based on the Rayleigh numbers assuming either basal heating ($\Delta T$) or internal heating (H).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. A2.** Flow regimes as a function of the unmolten and partially molten domains Rayleigh numbers ($Ra_{UM}$ and $Ra_{PM}$), after 20 Myr of basal heating of a 45 km thick continental crust containing white buoyant and low-viscosity material and black heavy and high-viscosity material (inclusions). The ambient domain remains motionless as long as $Ra_{PM} < 200$ (with at most local segregation of some inclusions one with respect to a neighbour), diapirism occurs when $200 < Ra_{PM} < 3000$, and the suspension regime initiates when $Ra_{UM} > 10$ (convection). The layering regime occurs at low $Ra_{UM}$ and $Ra_{PM} > 3000$: convection is not vigorous enough so that the buoyant inclusions can stack at ca. 20 km depth (after Louis-Napoléon et al., 2022).</figcaption></figure>