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

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.

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-migmatitic 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 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.

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):

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_p 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:

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

where K_eff is the consistency (in kg m⁻¹ s^−2+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_i is a material phase dependent pre-factor, with C_i the volume fraction of each phase i. The numerical term [$0.25×(0.75)−1n$] 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:

where α is the thermal expansion (α = 3 × 10⁻⁵ K⁻¹), T_ref is the temperature imposed at the top of the domain, and ρ_ref is a local phase-dependent density which depends on C_i (phase i volume fraction), and ρ_refi is a reference density for each material phase (see definition below).

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. The ambient medium, the white layers and the black layers are denoted hereafter as material “phases” which are identified by volume fractions C₁, C₂ and C₃, respectively.

The initial temperature linearly increases from T_c = 300 °C at −10 km depth to T_h = 600 °C at −45 km depth. At the onset of the numerical experiment, the basal temperature is switched to T_h+ = 1000 °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).

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¹⁹ Pa.s (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.

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 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 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).

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¹⁸ 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.

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.

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:

with H_T half of the crust’s height, ρ density, α the thermal expansion coefficient, ΔT_UM and ΔT_PM the temperature gradient across the unmolten (UM) and partially molten (PM) crustal domains, κ_UM and κ_PM the corresponding thermal conductivities, and K_eff^UP and κ_eff^PM the corresponding consistency (equivalent to the inverse of a viscosity depending on its power law exponent).