EDP Sciences logo
Submit your paper
Search
Menu
Home All issues Volume 195 (2024) BSGF - Earth Sci. Bull., 195 (2024) 21 Full HTML
Open Access
Issue
BSGF - Earth Sci. Bull.
Volume 195, 2024
Article Number 21
Number of page(s) 15
DOI https://doi.org/10.1051/bsgf/2024017
Published online 15 October 2024

© A. Louis-Napoléon et al., Published by EDP Sciences 2024

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 F1 and F2 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 F1/F2 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 F1/F2 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.

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

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

ρrefCp[∂T∂t+U⋅∇T]=−∇.(ρrefCp(T)κ(T)∇T)+Hr(1)

where T is temperature (K), ρ is the local density (kg m−3), U the local velocity (m.s−1), κ is the thermal diffusivity (m2 s−1), H is the radiogenic heat production (W m−3), and Cp is the heat capacity (J kg−1 K−1). 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:

ρref∂U∂t+ρrefU⋅∇U=−∇P+ρg+∇⋅[μ(∇U+(∇U)T)](2)

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

μ= Keff(T)ε˙1n−1 with Keff=∑i=13Ci [0.25×(0.75Ai)−1n×exp(QnRT)](3)

where Keff is the consistency (in kg m−1 s−2+1/n) and corresponds to a dynamic viscosity when n = 1. ε˙  has a minimal value set to 10−16 s−1, R = 8.314 J mol−1 K−1 is the universal gas constant, and Ai is a material phase dependent pre-factor, with Ci 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×105 J mol−1, 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/m3) and viscosity (Δµ ranging over 105 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:

ρ=ρref×[1−α(T−Tref)] with ρref=∑i=13Ciρrefi(4)

where α is the thermal expansion (α = 3 × 10−5 K−1), Tref is the temperature imposed at the top of the domain, and ρref is a local phase-dependent density which depends on Ci (phase i volume fraction), and ρ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. The ambient medium, the white layers and the black layers are denoted hereafter as material “phases” which are identified by volume fractions C1, C2 and C3, respectively.

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

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

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 ∼ 1019 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 ∼ 1019 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.

thumbnail 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 1016 to 1023 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 ∼ 1019 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
thumbnail 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.

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 1023 at ca. 650°C to 1016 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).

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

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

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.

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:

RaUM=2(HT2)2κUM(ραgΔTUM(HT2)2Κ~effUM)n,RaPM=ραgΔTPM(HT2)3κPMκ~effPM(5)

thumbnail 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 (ΔT) or internal heating (H).
thumbnail Fig. A2

Flow regimes as a function of the unmolten and partially molten domains Rayleigh numbers (RaUM and RaPM), 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 RaPM < 200 (with at most local segregation of some inclusions one with respect to a neighbour), diapirism occurs when 200 < RaPM <3000, and the suspension regime initiates when RaUM > 10 (convection). The layering regime occurs at low RaUM and RaPM > 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).

with HT half of the crust’s height, ρ density, α the thermal expansion coefficient, ΔTUM and ΔTPM the temperature gradient across the unmolten (UM) and partially molten (PM) crustal domains, κUM and κPM the corresponding thermal conductivities, and KeffUP and κeffPM the corresponding consistency (equivalent to the inverse of a viscosity depending on its power law exponent).

References

  • Arzi AA. 1978. Critical phenomena in the rheology of partially melted rocks. Tectonophysics 44: 173–184. [CrossRef] [Google Scholar]
  • Avigad D. 1998. High-pressure metamorphism and cooling on SE Naxos (Cyclades, Greece). Eur J Mineral 10: 1309–1320. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Brown M. 2001. Orogeny, migmatites and leucogranites: a review. J Earth Syst Sci 110: 313–336. [CrossRef] [Google Scholar]
  • 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 [CrossRef] [Google Scholar]
  • Brun J-P. 1983. L’origine des domes gneissiques; modèles et tests. Bull. Société Géologique Fr. S7-XXV: 219–228. [Google Scholar]
  • Brun J-P., Gapais D, Le Theoff B. 1981. The mantled gneiss domes of Kuopio (Finland): interfering diapirs. Tectonophysics 74: 283–304. [CrossRef] [Google Scholar]
  • Brun J-P., Sokoutis D, Van Den Driessche F J. 1994. Analogue modeling of detachment fault systems and core complexes. Geology 22: 319–322. [CrossRef] [Google Scholar]
  • Buck WR. 1991. Modes of continental lithospheric extension. J Geophys Res Solid Earth 96: 20161–20178. [CrossRef] [Google Scholar]
  • Buick IS. 1991. The late Alpine evolution of an extensional shear zone, Naxos, Greece. J Geol Soc 148: 93–103. [CrossRef] [Google Scholar]
  • Buick IS, Holland TJB. 1989. The P-T path associated with crustal extension, Naxos, Cyclades, Greece. Geol Soc Lond Spec Publ 43: 365–369. [CrossRef] [Google Scholar]
  • 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 [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Carslaw HS, Jaeger JC. 1959. Conduction of heat in solids, 2nd ed. Oxford: Clarendon Press. [Google Scholar]
  • Chen Y, Morgan WJ. 1990. A nonlinear rheology model for mid-ocean ridge axis topography. J Geophys Res 95: 17583. [CrossRef] [Google Scholar]
  • Collins WJ. 1989. Polydiapirism of the Archean Mount Edgar Batholith, Pilbara Block, Western Australia. Precambrian Res 43: 41–62. [CrossRef] [Google Scholar]
  • Collins WJ, Murphy JB, Blereau E, Huang H-Q. 2021. Water availability controls crustal melting temperatures. Lithos 402-403: 106351. [CrossRef] [Google Scholar]
  • Coney PJ, Harms TA. 1984. Cordilleran metamorphic core complexes: Cenozoic extensional relics of Mesozoic compression. Geology 12: 550. [CrossRef] [Google Scholar]
  • Cruden AR, Koyi H, Schmeling H. 1995. Diapiric basal entrainment of mafic into felsic magma. Earth Planet Sci Lett 131: 321–340. [CrossRef] [Google Scholar]
  • Davis GH. 1983. Shear-zone model for the origin of metamorphic core complexes. Geology 11: 342. [CrossRef] [Google Scholar]
  • 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 [CrossRef] [Google Scholar]
  • 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 [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Eskola PE. 1948. The problem of mantled gneiss domes. Q J Geol Soc 104: 461–476. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Hacker BR, Abers GA, Peacock SM. 2003. Subduction factory 1. Theoretical mineralogy, densities, seismic wave speeds, and H2O contents: subduction zone mineralogy and physical properties. J Geophys Res Solid Earth 108. https://doi.org/10.1029/2001JB001127 [Google Scholar]
  • 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 [CrossRef] [Google Scholar]
  • Jansen JB. 1973. Geological Map of Naxos. [Google Scholar]
  • Jansen JBH, Schuiling RD. 1976. Metamorphism on Naxos; petrology and geothermal gradients. Am J Sci 276: 1225–1253. [CrossRef] [Google Scholar]
  • Jolivet L, Brun J-P. 2010. Cenozoic geodynamic evolution of the Aegean. Int J Earth Sci 99: 109–138. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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 [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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 [Google Scholar]
  • 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. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Lister GS, Banga G, Feenstra A. 1984. Metamorphic core complexes of Cordilleran type in the Cyclades, Aegean Sea, Greece. Geology 12: 221. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [Google Scholar]
  • Martin D, Nokes R. 1989. A fluid-dynamical study of crystal settling in convecting magmas. J Petrol 30: 1471–1500. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Myers JS, Watkins KP. 1985. Origin of granite-greenstone patterns, Yilgarn Block, Western Australia. Geology 13: 778. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Peillod A, Majka J, Ring U, Drüppel K, Patten C, Karlsson A, Włodek A, Tehler 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Peillod A, Tehler 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 [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Rabillard A, Jolivet L, Arbaret L, Bessière E, Laurent V, Menant A, Augier R, Beaudoin A. 2018. Synextensional 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 [Google Scholar]
  • Ramberg H. 1981. The role of gravity in orogenic belts. Geol Soc Lond Spec Publ 9: 125–140. [CrossRef] [Google Scholar]
  • Ranalli G. 1995. Rheology of the Earth, Springer. ed. [Google Scholar]
  • 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. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Rey PF, Teyssier C, Kruckenberg SC, Whitney DL. 2011. Viscous collision in channel explains double domes in metamorphic core complexes. Geology 39: 387–390. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Roscoe R. 1952. The viscosity of suspensions of rigid spheres. Br J Appl Phys 3: 267–269. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Rosenberg CL, Handy MR. 2005. Experimental deformation of partially melted granite revisited: implications for the continental crust. J Metamorph Geol 23: 19–28. [CrossRef] [Google Scholar]
  • 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 [CrossRef] [Google Scholar]
  • Sawyer EW. 1994. Melt segregation in the continental crust. Geology 22: 1019–1022. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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 [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Talbot CJ. 1979. Infrastructural migmatitic upwelling in East Greenland interpreted as thermal convective structures. Precambrian Res 8: 77–93. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Turcotte D, Schubert G. 2014. Geodynamics:, 3rd ed. Cambridge University Press. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [Google Scholar]
  • van der Molen F I, Paterson MS. 1979. Experimental deformation of partially-melted granite. Contrib Mineral Petrol 70: 299–318. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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 [CrossRef] [Google Scholar]
  • 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 [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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 [Google Scholar]
  • 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 [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Vanderhaeghe O, Teyssier C. 2001. Partial melting and flow of orogens. Tectonophysics 342: 451–472. [CrossRef] [Google Scholar]
  • Vielzeuf D, Holloway JR. 1988. Experimental determination of the fluid-absent melting relations in the pelitic system. Contrib Mineral Petrol 98: 257–276. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • Weinberg RF. 1997. Diapir-driven crustal convection: decompression melting, renewal of the magma source and the origin of nested plutons. Tectonophysics 271: 217–229. [CrossRef] [Google Scholar]
  • Weinberg RF, Hasalová, P. 2015. Water-fluxed melting of the continental crust: a review. Lithos 212-215: 158–188. [CrossRef] [Google Scholar]
  • Weinberg RF, Schmeling H. 1992. Polydiapirs: multiwavelength gravity structures. J Struct Geol 14: 425–436. [CrossRef] [Google Scholar]
  • Whitney DL, Teyssier C, Vanderhaeghe O. 2004. Gneiss domes and crustal flow. Gneiss Domes Orogeny 380: 15. [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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 [Google Scholar]
  • 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. [CrossRef] [Google Scholar]
  • 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. [CrossRef] [Google Scholar]

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.

All Tables

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

In the text

All Figures

thumbnail 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).
In the text
thumbnail 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.
In the text
thumbnail 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 1016 to 1023 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 ∼ 1019 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
In the text
thumbnail 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.
In the text
thumbnail 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).
In the text
thumbnail 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 (ΔT) or internal heating (H).
In the text
thumbnail Fig. A2

Flow regimes as a function of the unmolten and partially molten domains Rayleigh numbers (RaUM and RaPM), 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 RaPM < 200 (with at most local segregation of some inclusions one with respect to a neighbour), diapirism occurs when 200 < RaPM <3000, and the suspension regime initiates when RaUM > 10 (convection). The layering regime occurs at low RaUM and RaPM > 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).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.