# Convection and segregation in heterogeneous orogenic crust with a VOF method – II: how to form migmatite domes

Aurélie Louis-Napoléon, Muriel Gerbault, Thomas Bonometti, Olivier

Vanderhaeghe, Roland Martin, Nathan Maury

<figure><img src="image_1.png" /></figure>

## To cite this version:

Aurélie Louis-Napoléon, Muriel Gerbault, Thomas Bonometti, Olivier Vanderhaeghe, Roland Martin, et al.. Convection and segregation in heterogeneous orogenic crust with a VOF method – II: how to form migmatite domes. Geophysical Journal International, 2024, 236 (1), pp.207-232. 10.1093/gji/ggad388. hal-04532187

HAL Id: hal-04532187

https://ut3-toulouseinp.hal.science/hal-04532187

Submitted on 4 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Convection and segregation in heterogeneous orogenic crust with a VOF method – II: how to form migmatite domes

Aurélie Louis-Napoléon<sup>1,2</sup> Muriel Gerbault<sup>1,2</sup> Thomas Bonometti,<sup>1</sup> Olivier Vanderhaeghe<sup>1,2</sup> Roland Martin<sup>2</sup> and Nathan Maury<sup>1</sup>

Accepted 2023 October 2. Received 2023 September 23; in original form 2022 November 29

## SUMMARY

Migmatites and granitic-gneisses exhumed in Archean to Phanerozoic segments are former partially molten crustal roots, display typical domes structures ranging in size from kilometres to decakilometres, and are often interpreted as resulting from the development of diapiric or convective gravitational instabilities. In previous work (part I), we determined various regimes of gravity-driven segregation, by considering a thick continental crust heated from below and containing melt related heterogeneities. These heterogeneities, represented by inclusions of distinct densities and viscosities with respect to the ambient partially molten material, can be entrained into convection cells (in the ‘suspension’ and ‘layering’ regimes) and/or accumulate as clusters (in the ‘layering’ and ‘diapirism’ regimes). Here we further investigate the specific conditions that allow for the formation and preservation of domes resulting from diapirism at the top of convective cells. We show that both the cessation of basal heating and the freezing of the buoyant inclusions density favour their stacking and preservation at ca. 15 km depth, within about 10 Myr. The buoyant inclusions form domes, 5–20 km in size, that also record several convective cycles at velocities ranging from 0.5–4 cm yr⁻<fnref n="1" />. 3-D models demonstrate their radial geometrical nature. The influence of the size and concentration of the inclusions is also assessed, complementing the characteristics of crustal heterogeneity in driving its differentiation and the formation of migmatite domes.

**Key words:** Composition and structure of the continental crust; Numerical modelling; Diapirism; Dynamics: gravity and tectonics; Heat generation and transport; Rheology: crust and lithosphere; Melt segregation; Migmatite domes.

## 1 INTRODUCTION

Domes cored by migmatites and granites, representing, respectively, former partially molten rocks and magmas, are ubiquitous structures of Archean to Phanerozoic crustal roots (Whitney et al. 2004; Clos et al. 2019). Many have been demonstrated to result from the development of gravitational instabilities (Platt 1980; Ramberg 1981b; Brun 1983; Van Kranendonk et al. 2004; Zibra et al. 2020). Structural analysis indicates that such domes develop from the crustal scale down to the decimetre scale with sometimes imbricated sub-domes into larger domes (Weinberg & Schmeling 1992; Burg & Vanderhaeghe 1993; Cruden et al. 1995; Vigneresse 1995; Ledru et al. 2001; Vidal et al. 2009; Brown 2010; Kruckenberg et al. 2011). In some cases, it has been argued that the structural, thermal and geochronological data records crustal-scale convection of the former partially molten orogenic roots (Collins et al. 1998; Van Kranendonk 2011; Riel et al. 2016; Vanderhaeghe et al. 2018).

In previous contributions, we presented how numerical modelling based on the VOF (volume-of-fluid) method was pertinent to capture the detailed structural evolution of a heterogeneous crustal system composed of fluids with strain-rate and temperature-dependent viscosities, heated from below (Louis-Napoléon et al. 2020, 2022). In these models, the heterogeneity of the orogenic crust was described by inclusions of several 100 m in size, having distinct density and viscosity with respect to the ambient bulk rock domain (called bulk equivalent domain ‘BED’). To justify this assumption, we used field observations suggesting that melt/crystal segregation in partially molten rocks occurs by percolation at the centimetre to decametre scale, as is attested by the presence of concordant to discordant granitic vein networks relative to the foliation of migmatites at this scale (Sawyer 1994; Brown et al. 1995; Vanderhaeghe & Teyssier 2001; Weinberg et al. 2015). These observations also corroborate theoretical studies predicting that felsic crustal melts can migrate via porosity waves at characteristic scales of a few hundred metres

<footnote n="1"><sup>1</sup>Institut de Mécanique des Fluides de Toulouse (IMFT), Université de Toulouse, CNRS, INPT, UPS, 31400 Toulouse, France. E-mail: aurelie.louisnapoleon@imft.fr</footnote>

<footnote n="2"><sup>2</sup>GET (UMR5563 CNRS, Univ. Paul Sabatier, IRD, CNES, OMP), 14 av. E. Belin, 31400 Toulouse, France</footnote>
# Convection and segregation in heterogeneous orogenic crust with a VOF method – II: how to form migmatite domes

Aurélie Louis-Napoléon<sup>1,2</sup> Muriel Gerbault<sup>1,2</sup> Thomas Bonometti,<sup>1</sup> Olivier Vanderhaeghe<sup>1,2</sup> Roland Martin<sup>2</sup> and Nathan Maury<sup>1</sup>

Accepted 2023 October 2. Received 2023 September 23; in original form 2022 November 29

## SUMMARY

Migmatites and granitic-gneisses exhumed in Archean to Phanerozoic segments are former partially molten crustal roots, display typical domes structures ranging in size from kilometres to decakilometres, and are often interpreted as resulting from the development of diapiric or convective gravitational instabilities. In previous work (part I), we determined various regimes of gravity-driven segregation, by considering a thick continental crust heated from below and containing melt related heterogeneities. These heterogeneities, represented by inclusions of distinct densities and viscosities with respect to the ambient partially molten material, can be entrained into convection cells (in the ‘suspension’ and ‘layering’ regimes) and/or accumulate as clusters (in the ‘layering’ and ‘diapirism’ regimes). Here we further investigate the specific conditions that allow for the formation and preservation of domes resulting from diapirism at the top of convective cells. We show that both the cessation of basal heating and the freezing of the buoyant inclusions density favour their stacking and preservation at ca. 15 km depth, within about 10 Myr. The buoyant inclusions form domes, 5–20 km in size, that also record several convective cycles at velocities ranging from 0.5–4 cm yr⁻<fnref n="1" />. 3-D models demonstrate their radial geometrical nature. The influence of the size and concentration of the inclusions is also assessed, complementing the characteristics of crustal heterogeneity in driving its differentiation and the formation of migmatite domes.

**Key words:** Composition and structure of the continental crust; Numerical modelling; Diapirism; Dynamics: gravity and tectonics; Heat generation and transport; Rheology: crust and lithosphere; Melt segregation; Migmatite domes.

## 1 INTRODUCTION

Domes cored by migmatites and granites, representing, respectively, former partially molten rocks and magmas, are ubiquitous structures of Archean to Phanerozoic crustal roots (Whitney et al. 2004; Clos et al. 2019). Many have been demonstrated to result from the development of gravitational instabilities (Platt 1980; Ramberg 1981b; Brun 1983; Van Kranendonk et al. 2004; Zibra et al. 2020). Structural analysis indicates that such domes develop from the crustal scale down to the decimetre scale with sometimes imbricated sub-domes into larger domes (Weinberg & Schmeling 1992; Burg & Vanderhaeghe 1993; Cruden et al. 1995; Vigneresse 1995; Ledru et al. 2001; Vidal et al. 2009; Brown 2010; Kruckenberg et al. 2011). In some cases, it has been argued that the structural, thermal and geochronological data records crustal-scale convection of the former partially molten orogenic roots (Collins et al. 1998; Van Kranendonk 2011; Riel et al. 2016; Vanderhaeghe et al. 2018).

In previous contributions, we presented how numerical modelling based on the VOF (volume-of-fluid) method was pertinent to capture the detailed structural evolution of a heterogeneous crustal system composed of fluids with strain-rate and temperature-dependent viscosities, heated from below (Louis-Napoléon et al. 2020, 2022). In these models, the heterogeneity of the orogenic crust was described by inclusions of several 100 m in size, having distinct density and viscosity with respect to the ambient bulk rock domain (called bulk equivalent domain ‘BED’). To justify this assumption, we used field observations suggesting that melt/crystal segregation in partially molten rocks occurs by percolation at the centimetre to decametre scale, as is attested by the presence of concordant to discordant granitic vein networks relative to the foliation of migmatites at this scale (Sawyer 1994; Brown et al. 1995; Vanderhaeghe & Teyssier 2001; Weinberg et al. 2015). These observations also corroborate theoretical studies predicting that felsic crustal melts can migrate via porosity waves at characteristic scales of a few hundred metres

<footnote n="1"><sup>1</sup>Institut de Mécanique des Fluides de Toulouse (IMFT), Université de Toulouse, CNRS, INPT, UPS, 31400 Toulouse, France. E-mail: aurelie.louisnapoleon@imft.fr</footnote>

<footnote n="2"><sup>2</sup>GET (UMR5563 CNRS, Univ. Paul Sabatier, IRD, CNES, OMP), 14 av. E. Belin, 31400 Toulouse, France</footnote>
(Connolly & Podladchikov 1998; Edmonds et al. 2019; Räss et al. 2019). The modelled inclusions aim at representing characteristic density and viscosity contrasts with respect to the 'BED'. When this orogenic crust was heated from below for about 20 Myr, the melting front propagated upwards and these inclusions became mobile, not only segregating individually but also being able to be entrained in global diapiric or convective flow. Our results led to the definition of four typical flow regimes that depend on the characteristic Rayleigh numbers of the unmolten ($Ra_{UM}$) and molten ($Ra_{PM}$) crustal domains (cf. details below). At low Rayleigh numbers, the inclusions are found to segregate individually due to their density contrasts. Diapirism initiates for a $Ra_{PM}$ exceeding $10^2$ and thermal convection occurs for a $Ra_{PM}$ over $10^3$. For intermediate values of $Ra_{PM}$ and $Ra_{UM}$, convection is associated with segregation of the inclusions that accumulate and form layers, whereas for high Rayleigh numbers, the inclusions are maintained in suspension in the convective cells. In the suspension regime, domes can form but they are reworked owing to convective motion, whereas in the layering regime, domes tend to flatten out as layers, with time.

In the present contribution, we further explore the conditions required for the preservation of dome structures in a convective context. We provide first a brief presentation of natural domes, describing the geodynamic settings within which they form and the smaller-scale melt transport processes they are associated to. The solved governing equations are then recalled, as well as our previous results (contribution part I), from which we select two specific configurations here. Then we present the results of new numerical experiments aimed at identifying additional conditions allowing for the preservation of domes within a convecting system. We test how the cessation of basal heating and the 'freezing' of the density of melt-related heterogeneities can actually favour dome formation. The parametric analysis is performed in 2-D, whereas the most pertinent configurations are extended in 3-D. The models insights on migmatite dome formation are then discussed, together with our approach's limitations and prospective work.

## 2 GEOLOGICAL CONSTRAINTS ON GNEISS DOMES AND MELT SEGREGATION-MIGRATION

### 2.1 On the geodynamic context for gneiss dome formation

Mantled-gneiss domes, first defined by Eskola (1948) to document the complex polyphased structures of nappes in the Scandinavian Shield, are ubiquitous features of exposed Archean to Phanerozoic crustal continental roots (Whitney et al. 2004; Clos et al. 2019, and references therein). Gneiss domes are typically cored by migmatites and granitoids, surrounded (i.e. mantled) by lower grade metamorphic rocks. They are defined by concentric trajectories of foliation or gneissic layering. Such domes are identified from kilometric to crustal scale and are characterized by a map-view axial ratio ranging between 1:1 and 3.5:1 irrespective of the size of the dome. Gneiss domes might be grouped in linear belts or distributed in domes and basins structures, in so-called granite-gneiss complexes of Archean, Proterozoic and Phanerozoic terranes (Fig. 1).

Mechanical models for the development of gneiss domes comprise (Burg et al. 2004; Yin 2004):

(i) Fold interferences caused by shortening in horizontal and orthogonal directions, eventually associated with thrust duplexes

(Ramsay 1967). These folds are characterized by mechanical coupling between the core gneisses and the surrounding metamorphic rocks. Such interferences are likely to occur in tectonically active, thickened crust coeval with active surface processes.

(ii) Detachment-controlled isostatic flow induced by vertical variations in viscosity (Buck 1991). It is characterized by a symptomatic superimposition of ductile and brittle structures, rapid exhumation associated with isothermal decompression of metamorphic rocks in the detachment's footwall, and by supradetachment basins in the hanging wall (Wernicke 1985).

(iii) Diapiric gravitational instabilities induced by an inverted density gradient (Ramberg 1981a; Weinberg & Schmeling 1992; Cruden et al. 1995). These are a favoured mechanism for gneiss dome formation when melt is present. Indeed, the density of silicate melt and magma is about 10 per cent lower than in typical continental crustal rocks (see discussion in Vanderhaeghe 2009). Diapirs are characterized by strain markers indicative of flattening at their top, dominant non-coaxial shearing along their flanks, coaxial radial stretching lineations or constriction in their core. Upwelling of the dome core is marked by cascading folds at its bordering boundaries (Kruckenberg et al. 2008; Clos et al. 2019). Diapiric flow is associated with near-isothermal decompression of the core rocks at a temperature exceeding the solidus, which might enhance dehydration-related partial melting (e.g. Duchêne et al. 2006).

These different models call upon distinct contributions of (i) surface processes related to plate tectonics induced topography and of (ii) buoyancy. Although petro-structural analyses provide characteristic length scales that are key to decipher the main forces at play, spacial arrangements point to the combined effects of gravity and tectonic forces (Vanderhaeghe et al. 2004; Kruckenberg et al. 2011). Foliation trajectories within migmatite terranes reveal second order domes, designated as subdomes. Spectacular examples comprise the Archean Pilbara gneiss domes (Collins 1989; Van Kranendonk et al. 2007; Clos et al. 2019; Roberts & Tikoff 2021), the Nassian Domain in Eburnean belt (Vidal et al. 2009) and the Vaasa migmatite complex of the Svecofennian belt (Chopin et al. 2020) for the Proterozoic period. Phanerozoic examples include the Carboniferous Velay dome of the Variscan belt (Ledru et al. 2001), the Jurassic Danba dome in the Songpan Ganze terrane (Billerot et al. 2017), the Eoocene Thor-Odin dome of the Shuswap Metamorphic Core Complex in the Canadian Cordillera (Vanderhaeghe et al. 1999) and the Miocene Naxos dome of the Aegean domain (Kruckenberg et al. 2011). The Naxos subdomes have been interpreted either as reflecting bidirectional isostatic flow (Rey et al. 2011) or as a combination of convective and diapiric gravitational instabilities (Vanderhaeghe et al. 2018).

Thermal convection has been further advocated based (i) on a low geothermal gradient identified on the exhumed root of the Ecuador arc (Riel et al. 2016) and (ii) on the cyclic geochronological record of zircon grains from the core of the Naxos migmatite dome with a period of 1–2 Myr (Vanderhaeghe et al. 2018). On the other hand, thermomechanical models using power-law crustal rheology over several tens of millions years have also shown that convection can reasonably develop within thickened hot orogenic crust (Cao et al. 2016), in contexts such as the present-day Altiplano Puna (Babeyko et al. 2002). This study deals with better identifying the thermomechanical conditions by which domes may form within a convecting orogenic crust.
(Connolly & Podladchikov 1998; Edmonds et al. 2019; Räss et al. 2019). The modelled inclusions aim at representing characteristic density and viscosity contrasts with respect to the 'BED'. When this orogenic crust was heated from below for about 20 Myr, the melting front propagated upwards and these inclusions became mobile, not only segregating individually but also being able to be entrained in global diapiric or convective flow. Our results led to the definition of four typical flow regimes that depend on the characteristic Rayleigh numbers of the unmolten ($Ra_{UM}$) and molten ($Ra_{PM}$) crustal domains (cf. details below). At low Rayleigh numbers, the inclusions are found to segregate individually due to their density contrasts. Diapirism initiates for a $Ra_{PM}$ exceeding $10^2$ and thermal convection occurs for a $Ra_{PM}$ over $10^3$. For intermediate values of $Ra_{PM}$ and $Ra_{UM}$, convection is associated with segregation of the inclusions that accumulate and form layers, whereas for high Rayleigh numbers, the inclusions are maintained in suspension in the convective cells. In the suspension regime, domes can form but they are reworked owing to convective motion, whereas in the layering regime, domes tend to flatten out as layers, with time.

In the present contribution, we further explore the conditions required for the preservation of dome structures in a convective context. We provide first a brief presentation of natural domes, describing the geodynamic settings within which they form and the smaller-scale melt transport processes they are associated to. The solved governing equations are then recalled, as well as our previous results (contribution part I), from which we select two specific configurations here. Then we present the results of new numerical experiments aimed at identifying additional conditions allowing for the preservation of domes within a convecting system. We test how the cessation of basal heating and the 'freezing' of the density of melt-related heterogeneities can actually favour dome formation. The parametric analysis is performed in 2-D, whereas the most pertinent configurations are extended in 3-D. The models insights on migmatite dome formation are then discussed, together with our approach's limitations and prospective work.

## 2 GEOLOGICAL CONSTRAINTS ON GNEISS DOMES AND MELT SEGREGATION-MIGRATION

### 2.1 On the geodynamic context for gneiss dome formation

Mantled-gneiss domes, first defined by Eskola (1948) to document the complex polyphased structures of nappes in the Scandinavian Shield, are ubiquitous features of exposed Archean to Phanerozoic crustal continental roots (Whitney et al. 2004; Clos et al. 2019, and references therein). Gneiss domes are typically cored by migmatites and granitoids, surrounded (i.e. mantled) by lower grade metamorphic rocks. They are defined by concentric trajectories of foliation or gneissic layering. Such domes are identified from kilometric to crustal scale and are characterized by a map-view axial ratio ranging between 1:1 and 3.5:1 irrespective of the size of the dome. Gneiss domes might be grouped in linear belts or distributed in domes and basins structures, in so-called granite-gneiss complexes of Archean, Proterozoic and Phanerozoic terranes (Fig. 1).

Mechanical models for the development of gneiss domes comprise (Burg et al. 2004; Yin 2004):

(i) Fold interferences caused by shortening in horizontal and orthogonal directions, eventually associated with thrust duplexes

(Ramsay 1967). These folds are characterized by mechanical coupling between the core gneisses and the surrounding metamorphic rocks. Such interferences are likely to occur in tectonically active, thickened crust coeval with active surface processes.

(ii) Detachment-controlled isostatic flow induced by vertical variations in viscosity (Buck 1991). It is characterized by a symptomatic superimposition of ductile and brittle structures, rapid exhumation associated with isothermal decompression of metamorphic rocks in the detachment's footwall, and by supradetachment basins in the hanging wall (Wernicke 1985).

(iii) Diapiric gravitational instabilities induced by an inverted density gradient (Ramberg 1981a; Weinberg & Schmeling 1992; Cruden et al. 1995). These are a favoured mechanism for gneiss dome formation when melt is present. Indeed, the density of silicate melt and magma is about 10 per cent lower than in typical continental crustal rocks (see discussion in Vanderhaeghe 2009). Diapirs are characterized by strain markers indicative of flattening at their top, dominant non-coaxial shearing along their flanks, coaxial radial stretching lineations or constriction in their core. Upwelling of the dome core is marked by cascading folds at its bordering boundaries (Kruckenberg et al. 2008; Clos et al. 2019). Diapiric flow is associated with near-isothermal decompression of the core rocks at a temperature exceeding the solidus, which might enhance dehydration-related partial melting (e.g. Duchêne et al. 2006).

These different models call upon distinct contributions of (i) surface processes related to plate tectonics induced topography and of (ii) buoyancy. Although petro-structural analyses provide characteristic length scales that are key to decipher the main forces at play, spacial arrangements point to the combined effects of gravity and tectonic forces (Vanderhaeghe et al. 2004; Kruckenberg et al. 2011). Foliation trajectories within migmatite terranes reveal second order domes, designated as subdomes. Spectacular examples comprise the Archean Pilbara gneiss domes (Collins 1989; Van Kranendonk et al. 2007; Clos et al. 2019; Roberts & Tikoff 2021), the Nassian Domain in Eburnean belt (Vidal et al. 2009) and the Vaasa migmatite complex of the Svecofennian belt (Chopin et al. 2020) for the Proterozoic period. Phanerozoic examples include the Carboniferous Velay dome of the Variscan belt (Ledru et al. 2001), the Jurassic Danba dome in the Songpan Ganze terrane (Billerot et al. 2017), the Eoocene Thor-Odin dome of the Shuswap Metamorphic Core Complex in the Canadian Cordillera (Vanderhaeghe et al. 1999) and the Miocene Naxos dome of the Aegean domain (Kruckenberg et al. 2011). The Naxos subdomes have been interpreted either as reflecting bidirectional isostatic flow (Rey et al. 2011) or as a combination of convective and diapiric gravitational instabilities (Vanderhaeghe et al. 2018).

Thermal convection has been further advocated based (i) on a low geothermal gradient identified on the exhumed root of the Ecuador arc (Riel et al. 2016) and (ii) on the cyclic geochronological record of zircon grains from the core of the Naxos migmatite dome with a period of 1–2 Myr (Vanderhaeghe et al. 2018). On the other hand, thermomechanical models using power-law crustal rheology over several tens of millions years have also shown that convection can reasonably develop within thickened hot orogenic crust (Cao et al. 2016), in contexts such as the present-day Altiplano Puna (Babeyko et al. 2002). This study deals with better identifying the thermomechanical conditions by which domes may form within a convecting orogenic crust.
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Examples of migmatite nested domes worldwide, with sizes ranging from 1 to 50 km: (a) the gneiss domes of the Archean Pilbara craton, Australia (after Collins 1989), (b) the internal structure of the Mt Edgar dome with geochronological data (after Roberts & Tikoff 2021; Vidal *et al.* 2009), (c) tonalite and granodiorite gneiss domes of the Palaeoproterozoic Nassian Domain, West Africa Craton, d) domes of the Vaasa migmatite complex in the Palaeoproterozoic Svecofennian belt (after Chopin *et al.* 2020), (e) the Velay dome of the Phanerozoic Variscan belt, (f) the Thor-Odin dome of the Shuswap Metamorphic Core Complex, Meso-Cenozoic Canadian Cordillera (after Vanderhaeghe *et al.* 1999) and (g) the Naxos dome of the Cenozoic Aegean domain (after Kruckenberg *et al.* 2011). Black lines with squares display low-angle detachment faults.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Examples of migmatite nested domes worldwide, with sizes ranging from 1 to 50 km: (a) the gneiss domes of the Archean Pilbara craton, Australia (after Collins 1989), (b) the internal structure of the Mt Edgar dome with geochronological data (after Roberts & Tikoff 2021; Vidal et al. 2009), (c) tonalite and granodiorite gneiss domes of the Palaeoproterozoic Nassian Domain, West Africa Craton, d) domes of the Vaasa migmatite complex in the Palaeoproterozoic Svecofennian belt (after Chopin et al. 2020), (e) the Velay dome of the Phanerozoic Variscan belt, (f) the Thor-Odin dome of the Shuswap Metamorphic Core Complex, Meso-Cenozoic Canadian Cordillera (after Vanderhaeghe et al. 1999) and (g) the Naxos dome of the Cenozoic Aegean domain (after Kruckenberg et al. 2011). Black lines with squares display low-angle detachment faults.</figcaption></figure>
## 2.2 Melt and magma transfer through the continental crust

Present day examples of large domains of partially molten continental crusts are documented by geophysical data, below the Central Andes and the Tibetan plateaus (e.g. Schilling et al. 2006; Caldwell et al. 2009). These magmas probably come from both the mantle and the crust in proportions that may vary depending on the geodynamic setting (e.g. Thompson & Connolly 1995). Furthermore, geophysical exploration methods at the crustal scale have been limited to resolutions greater than several kilometres, hence they cannot help evaluate the total amount of currently active proportions of partial melting. However, dense seismic array experiments such as those deployed recently at Long Beach/California (Bianco et al. 2019) or in the Pyrenees/France (Lehujeur & Chevrot 2020) were able to obtain high resolution tomographic images and detected the presence of fluids or gas down to approximately 10–15 km depth. Such teleseismic or local tomography seismic techniques aim at being applied to other areas and combined with gravity or magnetic data inversions (Chen et al. 2018; Spang et al. 2021), so that it should become possible in the mid-term to better estimate densities and fluid or melt fractions at mid-crustal depths.

A variety of mechanisms appear to be involved in the transfer of melt and magma through the crust, including (i) pervasive or localized porous flow by percolation through a solid mineral matrix (Leitch & Weinberg 2002; Bachmann & Huber 2016), (ii) mass transfer of magmas and/or partially molten rocks by diapiric flow (Ramberg 1981a; Schmeling et al. 2019) or (iii) localized transfer via dykes intruding the host rocks (Petford 1996). Some authors propose that melt transfer through the crust is achieved essentially by dyking (Annen & Sparks 2002), others favour percolation through magmatic columns (Bachmann & Huber 2016), but actually all these mechanisms might act in concert at different crustal levels (Vanderhaeghe 2009). Structural and petrological analyses of migmatites as well as analytical investigation of the compaction of a liquid/solid system, indicate that percolation is efficient at the grain scale and is associated with silicate melt segregation at the centimetre to decametre scales as granitic veins networks (Brown 1995; Sawyer et al. 1999; Vanderhaeghe & Teysnier 2001; Kisters et al. 2009; Weinberg et al. 2015). Melt migration out of the zone of partial melting might be impacted by heterogeneous deformation as indicated by the specific orientations of vein networks rooted in migmatites and feeding granitic plutons (Solar & Brown 2001; Vanderhaeghe 2009), which in turn display typical thicknesses of kilometric size (de Saint Blanquat et al. 2011). At shallower depths, crystallization of granitic magma is associated with preferred fractionation of the volatile phase within the last melt drops, and marks the transition from magmatic to hydrothermal fluid systems, which is an irreversible process (Balen & Broska 2011).

Crustal melting and melt transfer in orogenic crust has indeed been debated for decades and is still dividing the community. This debate may be sketched as two extreme views: one that sees crustal melting as a product that migrates buoyantly towards the surface via the porous crustal media, leaving a residual denser lower crust behind (e.g. Cashman et al. 2017), and one that considers that the entire partially molten mixture is mobilized, itself subject to gravity driven flow at a scale ranging from the grain size to the crustal breadth (Vanderhaeghe 2009). Our approach here rather stems from this latter view, as already discussed in part I (Louis-Napoléon et al. 2022): it ensues from previous theoretical studies of melt flow through porous partially molten continental crust, that obtain a typical scale of a few hundred metres for melt pockets resulting from

porosity waves propagation (Edmonds et al. 2019; Räss et al. 2019), a value which is also consistent with field analyses mentioned above (Brown 1995; Kisters et al. 2009; Brown 2010; Toé et al. 2013). This set of information drove our choice of defining inclusions of several hundred metres in size that would actually represent a mixture of partially molten material resulting from melt segregation, propagation and migration processes at the subscale. We chose a diameter size of 600 m, so as to approximate the fact that partial melting occurs without melt extraction beyond this scale. In order to depict the influence of lithological heterogeneities, we consider two end-members, namely inclusions that are either denser and more viscous than the ambient fluid, which mimic mafic lithologies, and inclusions that are lighter and less viscous than the ambient fluid, which mimic felsic lithologies. Hence, as in Louis-Napoléon et al. (2022), these inclusions characterize the heterogeneity scale of the partially molten crust.

# 3 PHYSICAL SETUP

## 3.1 Setup and qualitative description of the flow

The physical domain, used here to represent a standard orogenic crust, has horizontal dimensions of 50 km × 50 km and a height of 35 km (a domain depth range [−45, −10] km with origin set at the Earth surface), similarly to Louis-Napoléon et al. (2022). Fig. 2 presents the schematic chronological stages of the system:

(i) Initial stage: the temperature inside the domain is set to increase linearly from $T_C = 300$ °C at depth −10 km to $T_H = 600$ °C at depth −45 km. At time $t = 0$ Myr, the basal temperature is switched from $T_H = 600$ °C to $T_H = 1000$ °C. First, heat diffuses from the base of the domain across the bulk equivalent domain (BED, referring to material phase 1).

(ii) Heating stage: as melting progresses upwards from the model's base, convection initiates. Inclusions of contrasted density and viscosity relative to the BED are set in the domain (black and white inclusions, for heavy, more viscous material phase 3, and light, less viscous material phase 2, respectively). Prior to melting, their density and viscosity evolve similarly to those of the bulk domain, with a standard power-law strain-rate and temperature dependent rheology. Upon melting, and above a critical melt fraction $M_S$, the material becomes Newtonian, with a viscosity and density dependent on the melt fraction (as in Louis-Napoléon et al. 2022). In that case, we can consider that, at the first order, the location of the $M_S$ boundary represents the melting front: it separates the unmolten medium (UM, where $M < M_S$) from the partially molten medium (PM where $M > M_S$). The equations and the numerical setup are further described in sections below.

(iii) Cooling stage: the thermal boundary condition at the base is progressively cooled down from $T_H = 1000$ °C to $T_H = 600$ °C from time $t_1$ to time $t_2$. The typical duration of this process spans 5–10 Myr, according to the estimation of how fast can an underlying mantle lithosphere thermally recover from a process of delamination or slab break-off; if heat diffusion dominates in a lithospheric mantle about 20 km thick, diffusion scaling provides $\tau = H^2/\kappa \sim 4 \times 10^8/10^{-6} \sim 13$ Myr. During this cooling stage in our models (between $t_1$ and $t_2$), thermal convection ceases, diapirism may occur, and the system eventually returns to rest. Note that during the cooling process, the materials revert to their unmelted densities and viscosities, except for the models labelled 'freezing' where the inclusions densities adopt another trend (as explained in Section 3.3).
## 2.2 Melt and magma transfer through the continental crust

Present day examples of large domains of partially molten continental crusts are documented by geophysical data, below the Central Andes and the Tibetan plateaus (e.g. Schilling et al. 2006; Caldwell et al. 2009). These magmas probably come from both the mantle and the crust in proportions that may vary depending on the geodynamic setting (e.g. Thompson & Connolly 1995). Furthermore, geophysical exploration methods at the crustal scale have been limited to resolutions greater than several kilometres, hence they cannot help evaluate the total amount of currently active proportions of partial melting. However, dense seismic array experiments such as those deployed recently at Long Beach/California (Bianco et al. 2019) or in the Pyrenees/France (Lehujeur & Chevrot 2020) were able to obtain high resolution tomographic images and detected the presence of fluids or gas down to approximately 10–15 km depth. Such teleseismic or local tomography seismic techniques aim at being applied to other areas and combined with gravity or magnetic data inversions (Chen et al. 2018; Spang et al. 2021), so that it should become possible in the mid-term to better estimate densities and fluid or melt fractions at mid-crustal depths.

A variety of mechanisms appear to be involved in the transfer of melt and magma through the crust, including (i) pervasive or localized porous flow by percolation through a solid mineral matrix (Leitch & Weinberg 2002; Bachmann & Huber 2016), (ii) mass transfer of magmas and/or partially molten rocks by diapiric flow (Ramberg 1981a; Schmeling et al. 2019) or (iii) localized transfer via dykes intruding the host rocks (Petford 1996). Some authors propose that melt transfer through the crust is achieved essentially by dyking (Annen & Sparks 2002), others favour percolation through magmatic columns (Bachmann & Huber 2016), but actually all these mechanisms might act in concert at different crustal levels (Vanderhaeghe 2009). Structural and petrological analyses of migmatites as well as analytical investigation of the compaction of a liquid/solid system, indicate that percolation is efficient at the grain scale and is associated with silicate melt segregation at the centimetre to decametre scales as granitic veins networks (Brown 1995; Sawyer et al. 1999; Vanderhaeghe & Teyssier 2001; Kisters et al. 2009; Weinberg et al. 2015). Melt migration out of the zone of partial melting might be impacted by heterogeneous deformation as indicated by the specific orientations of vein networks rooted in migmatites and feeding granitic plutons (Solar & Brown 2001; Vanderhaeghe 2009), which in turn display typical thicknesses of kilometric size (de Saint Blanquat et al. 2011). At shallower depths, crystallization of granitic magma is associated with preferred fractionation of the volatile phase within the last melt drops, and marks the transition from magmatic to hydrothermal fluid systems, which is an irreversible process (Balen & Broska 2011).

Crustal melting and melt transfer in orogenic crust has indeed been debated for decades and is still dividing the community. This debate may be sketched as two extreme views: one that sees crustal melting as a product that migrates buoyantly towards the surface via the porous crustal media, leaving a residual denser lower crust behind (e.g. Cashman et al. 2017), and one that considers that the entire partially molten mixture is mobilized, itself subject to gravity driven flow at a scale ranging from the grain size to the crustal breadth (Vanderhaeghe 2009). Our approach here rather stems from this latter view, as already discussed in part I (Louis-Napoléon et al. 2022): it ensues from previous theoretical studies of melt flow through porous partially molten continental crust, that obtain a typical scale of a few hundred metres for melt pockets resulting from

porosity waves propagation (Edmonds et al. 2019; Räss et al. 2019), a value which is also consistent with field analyses mentioned above (Brown 1995; Kisters et al. 2009; Brown 2010; Toé et al. 2013). This set of information drove our choice of defining inclusions of several hundred metres in size that would actually represent a mixture of partially molten material resulting from melt segregation, propagation and migration processes at the subscale. We chose a diameter size of 600 m, so as to approximate the fact that partial melting occurs without melt extraction beyond this scale. In order to depict the influence of lithological heterogeneities, we consider two end-members, namely inclusions that are either denser and more viscous than the ambient fluid, which mimic mafic lithologies, and inclusions that are lighter and less viscous than the ambient fluid, which mimic felsic lithologies. Hence, as in Louis-Napoléon et al. (2022), these inclusions characterize the heterogeneity scale of the partially molten crust.

# 3 PHYSICAL SETUP

## 3.1 Setup and qualitative description of the flow

The physical domain, used here to represent a standard orogenic crust, has horizontal dimensions of 50 km × 50 km and a height of 35 km (a domain depth range [−45, −10] km with origin set at the Earth surface), similarly to Louis-Napoléon et al. (2022). Fig. 2 presents the schematic chronological stages of the system:

(i) Initial stage: the temperature inside the domain is set to increase linearly from $T_C = 300$ °C at depth −10 km to $T_H = 600$ °C at depth −45 km. At time $t = 0$ Myr, the basal temperature is switched from $T_H = 600$ °C to $T_H = 1000$ °C. First, heat diffuses from the base of the domain across the bulk equivalent domain (BED, referring to material phase 1).

(ii) Heating stage: as melting progresses upwards from the model's base, convection initiates. Inclusions of contrasted density and viscosity relative to the BED are set in the domain (black and white inclusions, for heavy, more viscous material phase 3, and light, less viscous material phase 2, respectively). Prior to melting, their density and viscosity evolve similarly to those of the bulk domain, with a standard power-law strain-rate and temperature dependent rheology. Upon melting, and above a critical melt fraction $M_S$, the material becomes Newtonian, with a viscosity and density dependent on the melt fraction (as in Louis-Napoléon et al. 2022). In that case, we can consider that, at the first order, the location of the $M_S$ boundary represents the melting front: it separates the unmolten medium (UM, where $M < M_S$) from the partially molten medium (PM where $M > M_S$). The equations and the numerical setup are further described in sections below.

(iii) Cooling stage: the thermal boundary condition at the base is progressively cooled down from $T_H = 1000$ °C to $T_H = 600$ °C from time $t_1$ to time $t_2$. The typical duration of this process spans 5–10 Myr, according to the estimation of how fast can an underlying mantle lithosphere thermally recover from a process of delamination or slab break-off; if heat diffusion dominates in a lithospheric mantle about 20 km thick, diffusion scaling provides $\tau = H^2/\kappa \sim 4 \times 10^8/10^{-6} \sim 13$ Myr. During this cooling stage in our models (between $t_1$ and $t_2$), thermal convection ceases, diapirism may occur, and the system eventually returns to rest. Note that during the cooling process, the materials revert to their unmelted densities and viscosities, except for the models labelled 'freezing' where the inclusions densities adopt another trend (as explained in Section 3.3).
<figure><img src="image_2.png" /><figcaption>**Figure 2.** (a) Scheme model setup. At time $t = 0$ Myr, the basal temperature is switched from $T_H = 600$ °C to $T_H + = 1000$ °C. Heat diffuses from the base across the bulk equivalent domain (BED, phase 1). As melting occurs, phase 1 properties and rheologies evolve (cf. Section 3.3) as two layers of unmolten (UM, blue domain, where the melt fraction $M \le M_s$) and partially molten material (PM, pink domain, where $M \ge M_s$). The latter contains heterogeneities—black and white deformable inclusions, corresponding to heavy more viscous material phase 3, and light less viscous material phase 2, respectively. Thermal convection may develop at this stage. Then, the basal temperature is progressively cooled down to $T_H = 600$ °C from time $t_1$ to time $t_2$. Convection dies out until the system eventually comes to rest. (b) Evolution of the local density and viscosity of each phase versus temperature (eqs 10, 5 and 8), at strain rate $\dot{\epsilon}_{min} = 10^{-16} s^{-1}$ for the UM viscosity. (c) Evolution of the local density in cases where we ‘freeze’ the light and heavy inclusions densities from time $t_1$ until the end of the simulation (cf. Section 3.3 and eq. 11).</figcaption></figure>

The heterogeneity of the partially molten domain is characterized by the presence of deformable inclusions of distinct densities and viscosities. These heterogeneous inclusions are supposed to have formed as a result from clustering melt percolation processes at the subscale, as ‘melt-pockets’ of typical sizes of ca. 600 m, a size that is justified in previous Section 2, and which will be discussed again in Section 5. They are labelled phases 2 and 3 (phase 1 corresponds to the surrounding BED).

## 3.2 System of equations

As in Louis-Napoléon et al. (2022), we consider the motion of a fluid with properties that may change due to partial melting, with the development of spatial heterogeneities once a certain amount of melt fraction is achieved. These heterogeneities are defined as distinct fluid phases $i$, identified by their volume fractions $C_i$ (index

$i = 1, 2, 3$ defines the phase and $\sum_{i=1}^3 C_i = 1$). More specifically, $C_i = 1$ inside phase $i$ and 0 outside, the three compositions never mix and we have sharp internal boundaries. The typical size of these internal boundaries is 1–2 grid cells. Viscosities and densities are thus a function of $C_i$ and are arithmetically averaged at the internal phase boundary. Hence the problem can be described by the following set of governing equations:

$$ \frac{\partial C_i}{\partial t} + \mathbf{U} \cdot \nabla C_i = -\nabla \cdot (\mathbf{U}_r C_r), \qquad (1) $$

$$ \nabla \cdot \mathbf{U} = 0, \qquad (2) $$

$$ \rho_{\text{ref}} \frac{\partial \mathbf{U}}{\partial t} + \rho_{\text{ref}} \mathbf{U} \cdot \nabla \mathbf{U} = -\nabla P + \bar{\rho} \mathbf{g} + \nabla \cdot [\hat{\mu}(\nabla \mathbf{U} + (\nabla \mathbf{U})^T)], \quad (3) $$

$$ \rho_{\text{ref}} C_p \left[ \frac{\partial T}{\partial t} + \mathbf{U} \cdot \nabla T \right] = \frac{1}{C^A} \left[ \nabla \cdot (\rho_{\text{ref}} C_p \kappa \nabla T) + H_r \right], \quad (4) $$
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Scheme model setup. At time $t = 0$ Myr, the basal temperature is switched from $T_H = 600$ °C to $T_H + = 1000$ °C. Heat diffuses from the base across the bulk equivalent domain (BED, phase 1). As melting occurs, phase 1 properties and rheologies evolve (cf. Section 3.3) as two layers of unmolten (UM, blue domain, where the melt fraction $M \le M_s$) and partially molten material (PM, pink domain, where $M \ge M_s$). The latter contains heterogeneities—black and white deformable inclusions, corresponding to heavy more viscous material phase 3, and light less viscous material phase 2, respectively. Thermal convection may develop at this stage. Then, the basal temperature is progressively cooled down to $T_H = 600$ °C from time $t_1$ to time $t_2$. Convection dies out until the system eventually comes to rest. (b) Evolution of the local density and viscosity of each phase versus temperature (eqs 10, 5 and 8), at strain rate $\dot{\epsilon}_{\min} = 10^{-16} s^{-1}$ for the UM viscosity. (c) Evolution of the local density in cases where we ‘freeze’ the light and heavy inclusions densities from time $t_1$ until the end of the simulation (cf. Section 3.3 and eq. 11).</figcaption></figure>

The heterogeneity of the partially molten domain is characterized by the presence of deformable inclusions of distinct densities and viscosities. These heterogeneous inclusions are supposed to have formed as a result from clustering melt percolation processes at the subscale, as ‘melt-pockets’ of typical sizes of ca. 600 m, a size that is justified in previous Section 2, and which will be discussed again in Section 5. They are labelled phases 2 and 3 (phase 1 corresponds to the surrounding BED).

## 3.2 System of equations

As in Louis-Napoléon et al. (2022), we consider the motion of a fluid with properties that may change due to partial melting, with the development of spatial heterogeneities once a certain amount of melt fraction is achieved. These heterogeneities are defined as distinct fluid phases $i$, identified by their volume fractions $C_i$ (index

$i = 1, 2, 3$ defines the phase and $\sum_{i=1}^3 C_i = 1$). More specifically, $C_i = 1$ inside phase $i$ and 0 outside, the three compositions never mix and we have sharp internal boundaries. The typical size of these internal boundaries is 1–2 grid cells. Viscosities and densities are thus a function of $C_i$ and are arithmetically averaged at the internal phase boundary. Hence the problem can be described by the following set of governing equations:

$$ \frac{\partial C_i}{\partial t} + \mathbf{U} \cdot \nabla C_i = -\nabla \cdot (\mathbf{U}_r C_r), \qquad (1) $$

$$ \nabla \cdot \mathbf{U} = 0, \qquad (2) $$

$$ \rho_{\text{ref}} \frac{\partial \mathbf{U}}{\partial t} + \rho_{\text{ref}} \mathbf{U} \cdot \nabla \mathbf{U} = -\nabla P + \bar{\rho} \mathbf{g} + \nabla \cdot [\hat{\mu}(\nabla \mathbf{U} + (\nabla \mathbf{U})^T)], \quad (3) $$

$$ \rho_{\text{ref}} C_p \left[ \frac{\partial T}{\partial t} + \mathbf{U} \cdot \nabla T \right] = \frac{1}{C^A} \left[ \nabla \cdot (\rho_{\text{ref}} C_p \kappa \nabla T) + H_r \right], \quad (4) $$
where **U**, **P**, **T** are the local velocity, pressure and temperature of the fluid, respectively, and **g** is the gravitational acceleration. The term $-\nabla(\mathbf{U}_r C_r)$ is artificially added to eq. (1) to reduce the effects of numerical smearing of the interface, and is specific to OpenFOAM's VOF method; $C_r = C_i \cdot (1 - C_i)$ for instance in a situation with only two phases, and $\mathbf{U}_r$, designated by Berberović et al. (2009) as a 'compression velocity', is evaluated at cell faces as a volume flux based on the maximum velocity magnitude in the interface region. This velocity is obtained from a face interpolation using the 'normalized variable diagram' (NVD approach Jakas et al. 1999).

In eq. (3), $\rho_{\text{ref}}$ and $\tilde{\rho}$ are the reference and the local density, respectively, whereas $\tilde{\mu}$ is the local dynamic viscosity. These quantities are detailed below. In eq. (4), $H_r = \sum_i C_i H_i^r$ is the total radioactive heating which depends on each phase $i$ and $H_i^r$ is the radioactive heating of phase $i$ (in W m⁻³), $\kappa(T)$ is the thermal diffusivity (in m² s⁻¹), and $C_p(T)$ is the heat capacity (in m² s⁻² K⁻¹). They are prescribed to depend on temperature according to Whittington et al. (2009), in a decoupled manner similar to (e.g. Gerya et al. 2008, cf. Appendix A2). We solve for latent heating via the term $C^A$ as well as for spatial variations in $C_p$ on the right-hand side of eq. (4), to the difference of part I (Louis-Napoléon et al. 2022). In particular, $C^A \approx 2$ when $T_{\text{sol}} \le T \le T_{\text{lig}}$ and $C^A = 1$ otherwise, as explained in Appendix A2.

### 3.3 Rheologies for the unmolten and partially molten domains

First, viscosities and densities are assigned to vary as a function of temperature and strain rate. As for the unmolten layer, the viscosity can be expressed as a consistency $K_{\text{eff}}$, which depends on a power-law relationship with temperature and strain rate following standard extrapolation from rock lab experiments (Chen & Morgan 1990):

$$ \tilde{\mu}^{\text{UM}} = K_{\text{eff}}^{\text{UM}}(T) [\max(\dot{\epsilon}_{\min}, \dot{\epsilon})]^{\frac{1}{n}-1}, \quad (5) $$

with

$$ K_{\text{eff}}^{\text{UM}}(T) = 0.25 \times 10^6 \times (0.75A)^{-\frac{1}{n}} \times \exp\left(\frac{Q}{nRT}\right), \quad (6) $$

where $T$ is the local temperature, $K_{\text{eff}}^{\text{UM}}$ is the consistency (in kg m⁻¹ s⁻² + 1/n) and corresponds to the dynamic viscosity when $n = 1$; $\dot{\epsilon} = \sqrt{(tr(\dot{\epsilon}^2) - tr^2(\dot{\epsilon})/3)/2}$, and $\dot{\epsilon} = (\nabla U + (\nabla U)^T)/2$ is the strain-rate tensor; $\dot{\epsilon}_{\min}$ is a minimal value for the strain rate set to $\dot{\epsilon}_{\min} = 10^{-16}$ s⁻¹; $R = 8.314$ J mol⁻¹ K⁻¹ is the universal gas constant; $Q = 1.54 \times 10^5$ J mol⁻¹ is the activation energy; $n = 2.3$ (parameters for quartz, Ranalli 1995). $A = \sum_{i=1}^3 C_i A_i$ is a constant pre-factor, where $A_i$ is defined for each phase $i$.

As of the rheology of partially molten crustal rocks, a number of studies have shown that it drastically changes from a given melt fraction threshold $M$ (Pinkerton & Stevenson 1992; Rosenberg & Handy 2005; Vanderhaeghe 2009). A few percent of melt are sufficient to cause a strength decrease of 2–3 orders of magnitude (Arzi 1978; Van der Molen & Paterson 1979; Rutter et al. 2006). On the other hand, the transition from magma to partially molten rocks is marked by a decrease in viscosity of several orders of magnitude (Roscoe 1952; Vetere et al. 2010; Schmeling et al. 2012). Furthermore, note that since there is no melt/solid extraction at the scale of our inclusions (but there might be at the subscale), their density is expected to decrease linearly with increasing melt fraction (see discussion in Vanderhaeghe 2009). These behaviours are formalized in the following way, with a melt fraction first defined for each phase as

$M^i = (T - T_{\text{sol}}^i)/(T_{\text{liq}}^i - T_{\text{sol}}^i)$. For simplicity then we choose equal solidus and liquidus temperatures for all phases, with $T_{\text{sol}} = 600$ °C and $T_{\text{liq}} = 1000$ °C, hence, as in Louis-Napoléon et al. (2022), we also have a single melt fraction $M$ such that:

$$ M = \begin{vmatrix} 0 & \text{if } T < T_{\text{sol}}, \\ \frac{T-T_{\text{sol}}}{T_{\text{liq}}-T_{\text{sol}}} & \text{if } T_{\text{sol}} \le T \le T_{\text{liq}}, \\ 1 & \text{if } T > T_{\text{liq}}. \end{vmatrix} \quad (7) $$

When $M$ exceeds $M_S$ (here $M_S = 0.3$), viscosities then switch to a Newtonian law which depends on $M$ and on a specific reference viscosity for each phase $\mu_i^0$, as proposed by (Pinkerton & Stevenson 1992, cf. Fig. 2b):

$$ \tilde{\mu}^{\text{PM}} = K_{\text{eff}}^{\text{PM}}(T) [\max(\dot{\epsilon}_{\min}, \dot{\epsilon})]^{\frac{1}{n}-1}, \quad (8) $$

with $n = 1$ and

$$ K_{\text{eff}}^{\text{PM}}(T) = \sum_{i=1}^{3} C_i \mu_i^0 \exp \left[ \left( 2.5 + \left( \frac{1-M}{M} \right)^{0.48} \right) (1-M) \right]. \quad (9) $$

For the densities in eq. (3), we use the Boussinesq approximation and we define an additional melt-dependent density (sic temperature-dependent, Fig. 2b) such as:

$$ \tilde{\rho} = \rho_{\text{ref}} \times [1 - \alpha(T - T_{\text{ref}})], \quad \text{and} \\ \rho_{\text{ref}} = \sum_{i=1}^{3} C_i [\rho_i^{S0}(1 - M) + \rho_i^{L0} M], \quad (10) $$

where $\alpha$ is the thermal expansion coefficient assumed constant ($\alpha = 3 \times 10^{-5}$ K⁻¹), $T_{\text{ref}}$ is the temperature imposed at the top of the domain (here $T_{\text{ref}} = T_C = 300$ °C at depth −10 km) and $\rho_{\text{ref}}$ is a local phase-dependent density which depends on $C_i$, the phase volume fraction, and $\rho_i^{S0}$, $\rho_i^{L0}$ the reference unmolten and molten densities for each fluid phase $i$, respectively.

In some tests below, we assume that during cooling the density difference between the light and heavy inclusions 'freezes', as a first order proxy to simulate the influence of crystallization processes. This assumption means that sufficient volatile extraction occurred at the subscale that cannot be recovered if reheating occurs. Indeed, magma crystallization is invoked to be associated with devolatilization and in particular the release of H₂O into the host rocks, which drives metasomatic alteration (Holloway 1976; Burnham 1979; Glazner 2019). This 'freezing' process is achieved by replacing $\rho_{\text{ref}}$ in eq. (10) by the following expression from time $t_1$ (when basal cooling is applied) until the end of the simulation:

$$ \rho_{\text{ref}} = C_1 \times [\rho_1^{S0}(1 - M) + \rho_1^{L0} M] \\ + \sum_{i=2}^{3} C_i [\rho_i^{S0}(1 - M^{\text{conv}}) + \rho_i^{L0} M^{\text{conv}}], \quad (11) $$

where $M^{\text{conv}}$ is the mean melt fraction in the convective region at time $t_1$. At this moment, the mean temperature in the convective region reaches about 800 °C, therefore $M^{\text{conv}} = (800 \text{ °C} - T_{\text{sol}})/(T_{\text{liq}} - T_{\text{sol}}) = 0.5$. The corresponding evolution of the inclusions and BED densities is displayed Fig. 2(c).

In summary, the unmelted viscosity of all materials is identical. Then upon melting they become Newtonian fluids where each material has a different viscosity, which varies with melt fraction. The same is true for densities. Note that upon cooling, materials can revert their viscosities and densities except for the models labelled 'freezing' for which the density difference between the light and heavy inclusions remains constant and becomes melt-independent. Finally, all materials have the same melting law, which means that at
where **U**, **P**, **T** are the local velocity, pressure and temperature of the fluid, respectively, and **g** is the gravitational acceleration. The term $-\nabla(\mathbf{U}_r C_r)$ is artificially added to eq. (1) to reduce the effects of numerical smearing of the interface, and is specific to OpenFOAM's VOF method; $C_r = C_i \cdot (1 - C_i)$ for instance in a situation with only two phases, and $\mathbf{U}_r$, designated by Berberović et al. (2009) as a 'compression velocity', is evaluated at cell faces as a volume flux based on the maximum velocity magnitude in the interface region. This velocity is obtained from a face interpolation using the 'normalized variable diagram' (NVD approach Jakas et al. 1999).

In eq. (3), $\rho_{\text{ref}}$ and $\tilde{\rho}$ are the reference and the local density, respectively, whereas $\tilde{\mu}$ is the local dynamic viscosity. These quantities are detailed below. In eq. (4), $H_r = \sum_i C_i H_i^r$ is the total radioactive heating which depends on each phase $i$ and $H_i^r$ is the radioactive heating of phase $i$ (in W m⁻³), $\kappa(T)$ is the thermal diffusivity (in m² s⁻¹), and $C_p(T)$ is the heat capacity (in m² s⁻² K⁻¹). They are prescribed to depend on temperature according to Whittington et al. (2009), in a decoupled manner similar to (e.g. Gerya et al. 2008, cf. Appendix A2). We solve for latent heating via the term $C^A$ as well as for spatial variations in $C_p$ on the right-hand side of eq. (4), to the difference of part I (Louis-Napoléon et al. 2022). In particular, $C^A \approx 2$ when $T_{\text{sol}} \le T \le T_{\text{lig}}$ and $C^A = 1$ otherwise, as explained in Appendix A2.

### 3.3 Rheologies for the unmolten and partially molten domains

First, viscosities and densities are assigned to vary as a function of temperature and strain rate. As for the unmolten layer, the viscosity can be expressed as a consistency $K_{\text{eff}}$, which depends on a power-law relationship with temperature and strain rate following standard extrapolation from rock lab experiments (Chen & Morgan 1990) :

$$ \tilde{\mu}^{\text{UM}} = K_{\text{eff}}^{\text{UM}}(T) [\max(\dot{\epsilon}_{\min}, \dot{\epsilon})]^{\frac{1}{n}-1}, \quad (5) $$

with

$$ K_{\text{eff}}^{\text{UM}}(T) = 0.25 \times 10^6 \times (0.75A)^{-\frac{1}{n}} \times \exp\left(\frac{Q}{nRT}\right), \quad (6) $$

where $T$ is the local temperature, $K_{\text{eff}}^{\text{UM}}$ is the consistency (in kg m⁻¹ s⁻² + 1/n) and corresponds to the dynamic viscosity when $n = 1$; $\dot{\epsilon} = \sqrt{(tr(\dot{\epsilon}^2) - tr^2(\dot{\epsilon})/3)/2}$, and $\dot{\epsilon} = (\nabla U + (\nabla U)^T)/2$ is the strain-rate tensor; $\dot{\epsilon}_{\min}$ is a minimal value for the strain rate set to $\dot{\epsilon}_{\min} = 10^{-16}$ s⁻¹; $R = 8.314$ J mol⁻¹ K⁻¹ is the universal gas constant; $Q = 1.54 \times 10^5$ J mol⁻¹ is the activation energy; $n = 2.3$ (parameters for quartz, Ranalli 1995). $A = \sum_{i=1}^3 C_i A_i$ is a constant pre-factor, where $A_i$ is defined for each phase $i$.

As of the rheology of partially molten crustal rocks, a number of studies have shown that it drastically changes from a given melt fraction threshold $M$ (Pinkerton & Stevenson 1992; Rosenberg & Handy 2005; Vanderhaeghe 2009). A few percent of melt are sufficient to cause a strength decrease of 2–3 orders of magnitude (Arzi 1978; Van der Molen & Paterson 1979; Rutter et al. 2006). On the other hand, the transition from magma to partially molten rocks is marked by a decrease in viscosity of several orders of magnitude (Roscoe 1952; Vetere et al. 2010; Schmeling et al. 2012). Furthermore, note that since there is no melt/solid extraction at the scale of our inclusions (but there might be at the subscale), their density is expected to decrease linearly with increasing melt fraction (see discussion in Vanderhaeghe 2009). These behaviours are formalized in the following way, with a melt fraction first defined for each phase as

$M^i = (T - T_{\text{sol}}^i)/(T_{\text{liq}}^i - T_{\text{sol}}^i)$. For simplicity then we choose equal solidus and liquidus temperatures for all phases, with $T_{\text{sol}} = 600$ °C and $T_{\text{liq}} = 1000$ °C, hence, as in Louis-Napoléon et al. (2022), we also have a single melt fraction $M$ such that:

$$ M = \begin{vmatrix} 0 & \text{if } T < T_{\text{sol}}, \\ \frac{T-T_{\text{sol}}}{T_{\text{liq}}-T_{\text{sol}}} & \text{if } T_{\text{sol}} \le T \le T_{\text{liq}}, \\ 1 & \text{if } T > T_{\text{liq}}. \end{vmatrix} \quad (7) $$

When $M$ exceeds $M_S$ (here $M_S = 0.3$), viscosities then switch to a Newtonian law which depends on $M$ and on a specific reference viscosity for each phase $\mu_i^0$, as proposed by (Pinkerton & Stevenson 1992, cf. Fig. 2b):

$$ \tilde{\mu}^{\text{PM}} = K_{\text{eff}}^{\text{PM}}(T) [\max(\dot{\epsilon}_{\min}, \dot{\epsilon})]^{\frac{1}{n}-1}, \quad (8) $$

with $n = 1$ and

$$ K_{\text{eff}}^{\text{PM}}(T) = \sum_{i=1}^{3} C_i \mu_i^0 \exp \left[ \left( 2.5 + \left( \frac{1-M}{M} \right)^{0.48} \right) (1-M) \right]. \quad (9) $$

For the densities in eq. (3), we use the Boussinesq approximation and we define an additional melt-dependent density (sic temperature-dependent, Fig. 2b) such as :

$$ \tilde{\rho} = \rho_{\text{ref}} \times [1 - \alpha(T - T_{\text{ref}})], \quad \text{and} \\ \rho_{\text{ref}} = \sum_{i=1}^{3} C_i [\rho_i^{S0}(1 - M) + \rho_i^{L0} M], \quad (10) $$

where $\alpha$ is the thermal expansion coefficient assumed constant ($\alpha = 3 \times 10^{-5}$ K⁻¹), $T_{\text{ref}}$ is the temperature imposed at the top of the domain (here $T_{\text{ref}} = T_C = 300$ °C at depth −10 km) and $\rho_{\text{ref}}$ is a local phase-dependent density which depends on $C_i$, the phase volume fraction, and $\rho_i^{S0}$, $\rho_i^{L0}$ the reference unmolten and molten densities for each fluid phase $i$, respectively.

In some tests below, we assume that during cooling the density difference between the light and heavy inclusions 'freezes', as a first order proxy to simulate the influence of crystallization processes. This assumption means that sufficient volatile extraction occurred at the subscale that cannot be recovered if reheating occurs. Indeed, magma crystallization is invoked to be associated with devolatilization and in particular the release of H₂O into the host rocks, which drives metasomatic alteration (Holloway 1976; Burnham 1979; Glazner 2019). This 'freezing' process is achieved by replacing $\rho_{\text{ref}}$ in eq. (10) by the following expression from time $t_1$ (when basal cooling is applied) until the end of the simulation:

$$ \rho_{\text{ref}} = C_1 \times [\rho_1^{S0}(1 - M) + \rho_1^{L0} M] + \sum_{i=2}^{3} C_i [\rho_i^{S0}(1 - M^{\text{conv}}) + \rho_i^{L0} M^{\text{conv}}], \quad (11) $$

where $M^{\text{conv}}$ is the mean melt fraction in the convective region at time $t_1$. At this moment, the mean temperature in the convective region reaches about 800 °C, therefore $M^{\text{conv}} = (800 \text{ °C} - T_{\text{sol}})/(T_{\text{liq}} - T_{\text{sol}}) = 0.5$. The corresponding evolution of the inclusions and BED densities is displayed Fig. 2(c).

In summary, the unmelted viscosity of all materials is identical. Then upon melting they become Newtonian fluids where each material has a different viscosity, which varies with melt fraction. The same is true for densities. Note that upon cooling, materials can revert their viscosities and densities except for the models labelled 'freezing' for which the density difference between the light and heavy inclusions remains constant and becomes melt-independent. Finally, all materials have the same melting law, which means that at
a given temperature, light and heavy inclusions and the BED have the same melt fraction but different densities and viscosities.

## 3.4 Numerical implementation

We use a VOF method implemented in the open source finite-volume code OpenFOAM, which we benchmarked in Louis-Napoléon et al. (2020) with Rayleigh-Taylor and Rayleigh-Bénard experiments from the literature. As in Louis-Napoléon et al. (2022), version 2106 of OpenFOAM was used, in which we developed our own solver multiMeltInterFoamv2 (with details in https://gitlab.com/AurelieLN/MultiMeltInterFoam). Coordinates x and z refer to the horizontal directions and y to the vertical direction. The mesh resolution used in the 2-D simulations is 250 × 175 for a domain size of 49.4 km × 35 km in the horizontal and vertical directions, respectively. In the 3-D simulations the mesh resolution is 250 × 175 × 250 with a regular grid spacing, and the domain size is 49.4 km × 35 km × 49.4 km in the (x, y, z)-directions, respectively.

The boundary conditions at the base and at the top of the model domain are the same in 2-D and in 3-D, with free slip at the base and no slip at the upper boundary ($y = -45$ km and $y = -10$ km, respectively), and we impose periodic lateral boundaries.

Initially, the inclusions occupy the entire spatial domain, with the same number of light and heavy inclusions (see Fig. 4). In two-dimensions, we use a staggered arrangement. The horizontal distance and vertical distance between the centres of two neighbouring inclusions are assumed equal and set to $L_t = 3r$. Here, $r = 600$ m. In 2-D, this leads to a theoretical volume fraction of all the inclusions (including light and heavy ones) equal to $\phi = \pi r^2 / L_t^2 \approx 0.35$. In three dimensions, we use a face-centred cubic lattice where the cube length is $L_t = (48)^{1/3}r \approx 3.63r$. Thus the volume fraction of all inclusions in 3-D is $\phi = 16\pi r^3 / (3L_t^3) \approx 0.35$, as in 2-D. Note that the effect of $\phi$ on the forecoming results is discussed in Appendix C.

The simulations ran on the regional cluster CALMIP https://www.calmip.univ-toulouse.fr/; for the 3-D simulations, they lasted 8 d on 36 cores for Configurations 1 and 2 (the configurations are described in the next section).

## 3.5 Results obtained in part I: four flow regimes

In Louis-Napoléon et al. (2022) we performed simulations that considered only a heating step from time 0 to about 20 Myr. Our parametric study led us to identify four main flow regimes (Fig. 3), found to occur depending on two Rayleigh numbers for the unmolten and partially molten domains, $Ra_{UM}$ and $Ra_{PM}$, respectively:

$$Ra_{UM} = \frac{2 \left(\frac{H_t}{2}\right)^2}{\kappa_{UM}} \left( \frac{\rho \alpha g \Delta T_{UM} \left(\frac{H_t}{2}\right)}{2 \bar{K}_{\text{eff}}^{\text{UM}}} \right)^n, \quad Ra_{PM} = \frac{\rho \alpha g \Delta T_{PM} \left(\frac{H_t}{2}\right)^3}{\kappa_{PM} \bar{K}_{\text{eff}}^{\text{PM}}}. \quad (12)$$

They characterize the thermal state of each domain, supposed to cover equally half the total crustal thickness, $H_t/2$ (hence here 22.5 km). Note that to each Ra correspond distinct temperature differences $\Delta T$, thermal diffusivities $\kappa$, and consistencies $K_{\text{eff}}$ (cf. details in Appendix A1).

The four flow regimes were described in part I. Let us briefly recall here the main difference between the suspension regime and the layering regime. Although convection occurs in both regimes, the light inclusions in the suspension regime cannot form a stagnant layer because the local velocity remains large in the region near the melting front where they could stack. This situation contrasts with that of the heavy and viscous inclusions, for which the bottom boundary is fixed and remains a low-velocity boundary layer. In this suspension regime, the actual low-velocity region (i.e. the boundary

layer) is located above the melting front. This result is coherent with a relatively small viscosity of the unmolten ambient rock (yielding a larger $Ra_{UM}$).

In contrast in the layering regime, the viscosity of the unmolten layer above the melting front is large enough to mimic a rigid lid: the viscosity gradient above the melting front is sharper, and the low-velocity boundary layer remains below the melting front. Thus, layering of the light inclusions at the top of the convective cells is favoured.

Note that only the layering regime allows for both convection and accumulation of the buoyant inclusions below the upper crust. They first stack as rounded shapes but then tend to flatten out over time and form a horizontal layer. Hence the purpose of this study part II to test complementary conditions that favour the formation and preservation of such domes.

# 4 RESULTS

We select now two distinct configurations (see parameters Table 1) within those displayed in Fig. 3: a Configuration 1 that lies in the ‘suspension regime’ where the light heterogeneities are entrained in the convective cells ($Ra_{UM} = 100$ and $Ra_{PM} = 3100$), and a Configuration 2 that lies in the ‘layering regime’ where the light heterogeneities can segregate above the convective cells (with $Ra_{UM} = 13$ and $Ra_{PM} = 3100$; cf. Fig. 3; Louis-Napoléon et al. 2022). Although there are several ways to switch from one regime to the other, thanks to Rayleigh’s dimensionalization properties, we control this switch here by modifying only the unmolten material viscosity. More specifically, the viscosity in Configuration 1 is set 2.3 times smaller than in Configuration 2 [since, according to eq. (5), the ratio of viscosities $\bar{\mu}^{\text{UM}}$ is $\propto K_{\text{eff}}^{\text{UM}}$ which is $\propto Ra_{UM}^{-1/n}$ via eq. (A1) with $n = 2.3$ here]. Then for each configuration, we proceed with a series of tests, displayed in Table 2. 2-D simulations are performed first: we test the onset and ending times of basal cooling and then the influence of freezing of the light inclusions upon cooling (details below). Finally, we select for each configuration the case considered as the most efficient in producing domes, and perform an equivalent 3-D simulation.

## 4.1 Configuration 1 – suspension regime

This configuration stands in the suspension regime according to Fig. 3, and has density and viscosity parameters displayed in Table 1. The density and viscosity of the inclusions (phases 2 and 3) get activated when $M > M_S$ (here $M_S = 0.3$, eqs 9, 10). Note that the effect of the specific value of $M_S$ was discussed in Louis-Napoléon et al. (2022; cf. their appendix A1).

### 4.1.1 Configuration 1 – No cooling

Here, basal heating is maintained over the model’s time duration, over 20 Myr (model $Config_1^\infty$). Snapshots in time are displayed Fig. 4(a). Thermal convection initiates around 4 Myr, and its vigour increases over time, with the heavy and light inclusions being continuously dragged into the convection cells. Although a majority of the heavy inclusions deposit at the base of the model domain, the light inclusions cannot segregate upwards. This vigorous convection regime generates large homogenized cells of characteristic width of about 20 km, function of the system’s thickness, heating time duration, and density and viscosity properties (cf. Rayleigh numbers determined in Louis-Napoléon et al. 2022).
a given temperature, light and heavy inclusions and the BED have the same melt fraction but different densities and viscosities.

## 3.4 Numerical implementation

We use a VOF method implemented in the open source finite-volume code OpenFOAM, which we benchmarked in Louis-Napoléon et al. (2020) with Rayleigh-Taylor and Rayleigh-Bénard experiments from the literature. As in Louis-Napoléon et al. (2022), version 2106 of OpenFOAM was used, in which we developed our own solver multiMeltInterFoamv2 (with details in https://gitlab.com/AurelieLN/MultiMeltInterFoam). Coordinates x and z refer to the horizontal directions and y to the vertical direction. The mesh resolution used in the 2-D simulations is 250 × 175 for a domain size of 49.4 km × 35 km in the horizontal and vertical directions, respectively. In the 3-D simulations the mesh resolution is 250 × 175 × 250 with a regular grid spacing, and the domain size is 49.4 km × 35 km × 49.4 km in the (x, y, z)-directions, respectively.

The boundary conditions at the base and at the top of the model domain are the same in 2-D and in 3-D, with free slip at the base and no slip at the upper boundary ($y = -45$ km and $y = -10$ km, respectively), and we impose periodic lateral boundaries.

Initially, the inclusions occupy the entire spatial domain, with the same number of light and heavy inclusions (see Fig. 4). In two-dimensions, we use a staggered arrangement. The horizontal distance and vertical distance between the centres of two neighbouring inclusions are assumed equal and set to $L_t = 3r$. Here, $r = 600$ m. In 2-D, this leads to a theoretical volume fraction of all the inclusions (including light and heavy ones) equal to $\phi = \pi r^2 / L_t^2 \approx 0.35$. In three dimensions, we use a face-centred cubic lattice where the cube length is $L_t = (48)^{1/3}r \approx 3.63r$. Thus the volume fraction of all inclusions in 3-D is $\phi = 16\pi r^3 / (3L_t^3) \approx 0.35$, as in 2-D. Note that the effect of $\phi$ on the forecoming results is discussed in Appendix C.

The simulations ran on the regional cluster CALMIP https://www.calmip.univ-toulouse.fr/; for the 3-D simulations, they lasted 8 d on 36 cores for Configurations 1 and 2 (the configurations are described in the next section).

## 3.5 Results obtained in part I: four flow regimes

In Louis-Napoléon et al. (2022) we performed simulations that considered only a heating step from time 0 to about 20 Myr. Our parametric study led us to identify four main flow regimes (Fig. 3), found to occur depending on two Rayleigh numbers for the unmolten and partially molten domains, $Ra_{UM}$ and $Ra_{PM}$, respectively:

$$Ra_{UM} = \frac{2 \left(\frac{H_t}{2}\right)^2}{\kappa_{UM}} \left( \frac{\rho \alpha g \Delta T_{UM} \left(\frac{H_t}{2}\right)}{2 \bar{K}_{\text{eff}}^{\text{UM}}} \right)^n, \quad Ra_{PM} = \frac{\rho \alpha g \Delta T_{PM} \left(\frac{H_t}{2}\right)^3}{\kappa_{PM} \bar{K}_{\text{eff}}^{\text{PM}}}. \quad (12)$$

They characterize the thermal state of each domain, supposed to cover equally half the total crustal thickness, $H_t/2$ (hence here 22.5 km). Note that to each Ra correspond distinct temperature differences $\Delta T$, thermal diffusivities $\kappa$, and consistencies $K_{\text{eff}}$ (cf. details in Appendix A1).

The four flow regimes were described in part I. Let us briefly recall here the main difference between the suspension regime and the layering regime. Although convection occurs in both regimes, the light inclusions in the suspension regime cannot form a stagnant layer because the local velocity remains large in the region near the melting front where they could stack. This situation contrasts with that of the heavy and viscous inclusions, for which the bottom boundary is fixed and remains a low-velocity boundary layer. In this suspension regime, the actual low-velocity region (i.e. the boundary

layer) is located above the melting front. This result is coherent with a relatively small viscosity of the unmolten ambient rock (yielding a larger $Ra_{UM}$).

In contrast in the layering regime, the viscosity of the unmolten layer above the melting front is large enough to mimic a rigid lid: the viscosity gradient above the melting front is sharper, and the low-velocity boundary layer remains below the melting front. Thus, layering of the light inclusions at the top of the convective cells is favoured.

Note that only the layering regime allows for both convection and accumulation of the buoyant inclusions below the upper crust. They first stack as rounded shapes but then tend to flatten out over time and form a horizontal layer. Hence the purpose of this study part II to test complementary conditions that favour the formation and preservation of such domes.

# 4 RESULTS

We select now two distinct configurations (see parameters Table 1) within those displayed in Fig. 3: a Configuration 1 that lies in the ‘suspension regime’ where the light heterogeneities are entrained in the convective cells ($Ra_{UM} = 100$ and $Ra_{PM} = 3100$), and a Configuration 2 that lies in the ‘layering regime’ where the light heterogeneities can segregate above the convective cells (with $Ra_{UM} = 13$ and $Ra_{PM} = 3100$; cf. Fig. 3; Louis-Napoléon et al. 2022). Although there are several ways to switch from one regime to the other, thanks to Rayleigh’s dimensionalization properties, we control this switch here by modifying only the unmolten material viscosity. More specifically, the viscosity in Configuration 1 is set 2.3 times smaller than in Configuration 2 [since, according to eq. (5), the ratio of viscosities $\bar{\mu}^{\text{UM}}$ is $\propto K_{\text{eff}}^{\text{UM}}$ which is $\propto Ra_{UM}^{-1/n}$ via eq. (A1) with $n = 2.3$ here]. Then for each configuration, we proceed with a series of tests, displayed in Table 2. 2-D simulations are performed first: we test the onset and ending times of basal cooling and then the influence of freezing of the light inclusions upon cooling (details below). Finally, we select for each configuration the case considered as the most efficient in producing domes, and perform an equivalent 3-D simulation.

## 4.1 Configuration 1 – suspension regime

This configuration stands in the suspension regime according to Fig. 3, and has density and viscosity parameters displayed in Table 1. The density and viscosity of the inclusions (phases 2 and 3) get activated when $M > M_S$ (here $M_S = 0.3$, eqs 9, 10). Note that the effect of the specific value of $M_S$ was discussed in Louis-Napoléon et al. (2022; cf. their appendix A1).

### 4.1.1 Configuration 1 – No cooling

Here, basal heating is maintained over the model’s time duration, over 20 Myr (model $Config_1^\infty$). Snapshots in time are displayed Fig. 4(a). Thermal convection initiates around 4 Myr, and its vigour increases over time, with the heavy and light inclusions being continuously dragged into the convection cells. Although a majority of the heavy inclusions deposit at the base of the model domain, the light inclusions cannot segregate upwards. This vigorous convection regime generates large homogenized cells of characteristic width of about 20 km, function of the system’s thickness, heating time duration, and density and viscosity properties (cf. Rayleigh numbers determined in Louis-Napoléon et al. 2022).
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Flow regimes obtained in Louis-Napoléon *et al.* (2022) as a function of the unmolten and partially molten domains Rayleigh numbers ($Ra_{UM}$ and $Ra_{PM}$), after 20 Myr of basal heating. The BED (bulk domain, phase 1) remains motionless as long as $Ra_{PM} < 200$, diapirism occurs when $200 < Ra_{PM} < 3000$, and the suspension regime initiates when $Ra_{UM} > 10$. The layering regime occurs at low $Ra_{UM}$ and $Ra_{PM} > 3000$: convection is not vigorous and the inclusions can stack at $ca. 20$ km depth. Here we explore the conditions for which dome formation is favoured in the two latter regimes noted Configurations 1 and 2.</figcaption></figure>

### 4.1.2 Configuration 1 – With cooling

Now, we progressively reduce the basal heating, starting from time $t_1 = 4.4$ Myr until time $t_2 = 9.4$ Myr (Fig. 4b) or 14.4 Myr (Fig. 4c) when the basal temperature returns to $T = 600$ °C. From time $t_1$ onwards, the light inclusions above the melting front are still being stretched and attracted into the convective flow, but some of them manage to accumulate below the rigid upper crust. Note, however, that no dome is observed in these cases. The heavy inclusions, in turn, keep settling down at the base of the physical domain and form a layer. This layer is more heterogeneous and less compact than in the case without cooling (Figs 4b and c versus Fig. 4a). Indeed, the local viscosity increases in the BED when cooling is applied and thus some interstitial BED fluid remains in between the heavy inclusions.

Complementary tests allow us to confirm the effect of cooling time with respect to the system's dynamics: applying it sooner (before 4.4 Myr) impedes the development of thermal convection, whereas applying it over a longer period delays the interruption of convection, hence delays the stacking of the light inclusions below the upper crust (cf. Figs 4b and c).

### 4.1.3 Configuration 1 – With cooling and freezing of the inclusions densities

Now, we impose that the density of the inclusions becomes melt-independent from time $t_1$. This option reflects the irreversibility of several crystallization processes (e.g. the extraction of most volatile fluids, cf. details in Section 2). Results are displayed in Fig. 4(d).

From time 9.4 Myr here, the buoyant inclusions start to form rounded shape clusters that accumulate below the rigid upper crust, around 15–20 km depth. These clusters display noticeable ribbon structures, that store more and more ‘white material’ over time (until 20 Myr). Their characteristic size is about 5 km.

The averaged vertical velocities displayed in Fig. 4(e) show oscillations inherent to cyclic convective instabilities, and the slowing down of the systems dynamics upon cooling the base of the crust from 4.4 Myr onwards. Although the ‘freezing densities’ model also includes basal cooling, its velocities decrease more than twice as fast as the same model without ‘freezing’ (the dashed green line compares with the dotted red line). We interpret this as the freezing of the inclusions densities favouring their accumulation
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Flow regimes obtained in Louis-Napoléon *et al.* (2022) as a function of the unmolten and partially molten domains Rayleigh numbers ($Ra_{UM}$ and $Ra_{PM}$), after 20 Myr of basal heating. The BED (bulk domain, phase 1) remains motionless as long as $Ra_{PM} < 200$, diapirism occurs when $200 < Ra_{PM} < 3000$, and the suspension regime initiates when $Ra_{UM} > 10$. The layering regime occurs at low $Ra_{UM}$ and $Ra_{PM} > 3000$: convection is not vigorous and the inclusions can stack at $ca. 20$ km depth. Here we explore the conditions for which dome formation is favoured in the two latter regimes noted Configurations 1 and 2.</figcaption></figure>

### 4.1.2 Configuration 1 – With cooling

Now, we progressively reduce the basal heating, starting from time $t_1 = 4.4$ Myr until time $t_2 = 9.4$ Myr (Fig. 4b) or 14.4 Myr (Fig. 4c) when the basal temperature returns to $T = 600$ °C. From time $t_1$ onwards, the light inclusions above the melting front are still being stretched and attracted into the convective flow, but some of them manage to accumulate below the rigid upper crust. Note, however, that no dome is observed in these cases. The heavy inclusions, in turn, keep settling down at the base of the physical domain and form a layer. This layer is more heterogeneous and less compact than in the case without cooling (Figs 4b and c versus Fig. 4a). Indeed, the local viscosity increases in the BED when cooling is applied and thus some interstitial BED fluid remains in between the heavy inclusions.

Complementary tests allow us to confirm the effect of cooling time with respect to the system’s dynamics: applying it sooner (before 4.4 Myr) impedes the development of thermal convection, whereas applying it over a longer period delays the interruption of convection, hence delays the stacking of the light inclusions below the upper crust (cf. Figs 4b and c).

### 4.1.3 Configuration 1 – With cooling and freezing of the inclusions densities

Now, we impose that the density of the inclusions becomes melt-independent from time $t_1$. This option reflects the irreversibility of several crystallization processes (e.g. the extraction of most volatile fluids, cf. details in Section 2). Results are displayed in Fig. 4(d).

From time 9.4 Myr here, the buoyant inclusions start to form rounded shape clusters that accumulate below the rigid upper crust, around 15–20 km depth. These clusters display noticeable ribbon structures, that store more and more ‘white material’ over time (until 20 Myr). Their characteristic size is about 5 km.

The averaged vertical velocities displayed in Fig. 4(e) show oscillations inherent to cyclic convective instabilities, and the slowing down of the systems dynamics upon cooling the base of the crust from 4.4 Myr onwards. Although the ‘freezing densities’ model also includes basal cooling, its velocities decrease more than twice as fast as the same model without ‘freezing’ (the dashed green line compares with the dotted red line). We interpret this as the freezing of the inclusions densities favouring their accumulation
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Configuration 1: Summary of the 2-D results. From top to bottom: (a) *Config1*∞, (b) *Config1*5−10, (c) *Config1*5−15 and (d) *Config1*F5−10 (see Tables 1 and 2 for parameters). For each case, snapshots display temperature isocontours in Celsius degrees and, on the left half of the domain the viscosity distribution, and on the right half, the light and heavy inclusions in white and black, respectively. (e) Time evolution of |*U*<sub>y</sub>| which denotes the spatial average of the absolute value of the vertical velocity.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Configuration 1: Summary of the 2-D results. From top to bottom: (a) *Config1*∞, (b) *Config1*5−10, (c) *Config1*5−15 and (d) *Config1*F5−10 (see Tables 1 and 2 for parameters). For each case, snapshots display temperature isocontours in Celsius degrees and, on the left half of the domain the viscosity distribution, and on the right half, the light and heavy inclusions in white and black, respectively. (e) Time evolution of |*U*<sub>y</sub>| which denotes the spatial average of the absolute value of the vertical velocity.</figcaption></figure>
<table><caption>**Table 1.** Parameters for Configurations 1 and 2: *i* is phase index ($i=1$ for the BED, $i=2$ for the light and less viscous inclusions and $i=3$ for the heavy more viscous inclusions); PM, partially molten; UM, unmolten material; $\rho_i^{L0}$ and $\rho_i^{S0}$ are the reference densities used in eq. (10); $H_i^r$ is radioactive heat production, $\mu_i^0$ and $A_i$ are viscosity pre-factors in eqs (9) and (6), respectively (for Configurations 1 or 2, labelled in exponent). Furthermore, $n=2.3$, $Q=1.54 \times 10^5$ J mol$^{-1}$, $R=8.314$ J mol$^{-1}$ K$^{-1}$, $\alpha=3 \times 10^{-5}$ K$^{-1}$, so that $Ra_{UM}=100$ and $Ra_{PM}=3100$ in Configuration 1 and $Ra_{UM}=13$ and $Ra_{PM}=3100$ in Configuration 2.</caption><thead><tr><th>Fluid</th><th>Density<br>(PM layer)</th><th>Density<br>(UM layer)</th><th>Production</th><th>Viscosity<br>(PM layer)</th><th>Viscosity pre-factor<br>(UM layer)</th><th>Viscosity pre-factor<br>(UM layer)</th></tr><tr><th><i>i</i></th><th><i>ρ</i><sub><i>i</i></sub><sup>L0</sup><br>(kg m<sup>-3</sup>)</th><th><i>ρ</i><sub><i>i</i></sub><sup>S0</sup><br>(kg m<sup>-3</sup>)</th><th><i>H</i><sub><i>i</i></sub><sup>r</sup><br>(μW m<sup>-3</sup>)</th><th><i>μ</i><sub><i>i</i></sub><sup>0</sup><br>(Pa s)</th><th>Config<sub>1</sub>: <i>A</i><sub><i>i</i></sub><sup>1</sup><br>(Pa<sup>-2.3</sup> s<sup>-1</sup>)</th><th>Config<sub>2</sub>: <i>A</i><sub><i>i</i></sub><sup>2</sup><br>(Pa<sup>-2.3</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1</td><td>2700</td><td>2800</td><td>1.75</td><td>8.7 × 10<sup>16</sup></td><td>3.2 × 10<sup>-4</sup></td><td>4.2 × 10<sup>-5</sup></td></tr><tr><td>2</td><td>2500</td><td>2800</td><td>10</td><td>8.7 × 10<sup>15</sup></td><td>3.2 × 10<sup>-4</sup></td><td>4.2 × 10<sup>-5</sup></td></tr><tr><td>3</td><td>2900</td><td>2800</td><td>0.9</td><td>8.7 × 10<sup>17</sup></td><td>3.2 × 10<sup>-4</sup></td><td>4.2 × 10<sup>-5</sup></td></tr></tbody></table>

<table><caption>**Table 2.** Configurations names (*Config*₁ and *Config*₂) distinguish options: the exponents indicate the onset and ending times of cooling at the bottom boundary ($\infty$ = no cooling); ‘F’ refers to freezing of the inclusions density from cooling time $t_1$; (3-D) indicates a 3-D simulation. The models produce convection characterized by a ‘Cycle time’, for example the convection period, and ‘Velocity’, for example the mean markers velocity range. The two last columns display convective cells temperature range (prior to the onset of cooling and measured from a non-exhaustive selection of markers), and the final size of the domes when the system has come to rest.</caption><thead><tr><th></th><th>Simulation name</th><th>Cooling period<br>[t<sub>1</sub>, t<sub>2</sub>] (in Myr)</th><th>Cycle time (τ)</th><th>Convection properties<br>Velocity (v<sub>rms</sub>)<br>[T<sub>min</sub>, T<sub>max</sub>; ΔT]</th><th>Domes<br>D<sub>20Myr</sub></th></tr></thead><tbody><tr><td rowspan="5">Configuration 1</td><td>Config<sub>1</sub><sup>∞</sup></td><td>[∞, ∞]</td><td>1.8 Myr</td><td>2–5 cm yr<sup>-1</sup></td><td>Ø;</td></tr><tr><td>Config<sub>1</sub><sup>5–10</sup></td><td>[4.4,9.4]</td><td>”</td><td>”</td><td>Ø;</td></tr><tr><td>Config<sub>1</sub><sup>5–15</sup></td><td>[4.4,14.4]</td><td>”</td><td>”</td><td>Ø;</td></tr><tr><td>Config<sub>1</sub><sup>F<sup>5</sup>–10</sup></td><td>[4.4,9.4]</td><td>”</td><td>”</td><td>10 km</td></tr><tr><td>Config<sub>1</sub><sup>F<sup>8</sup>–13</sup>(3-D)</td><td>[8.0,13.0]</td><td>2 Myr</td><td>2–3 cm yr<sup>-1</sup></td><td>[700, 900; 200] °C</td><td>10 km</td></tr><tr><td rowspan="5">Configuration 2</td><td>Config<sub>2</sub><sup>∞</sup></td><td>[∞, ∞]</td><td>2 Myr</td><td>0.5–1 cm yr<sup>-1</sup></td><td>[700, 850; 150] °C</td><td>5–10 km</td></tr><tr><td>Config<sub>2</sub><sup>7–12</sup></td><td>[7.0, 12.0]</td><td>”</td><td>”</td><td>”</td><td>5–10 km</td></tr><tr><td>Config<sub>2</sub><sup>F<sup>7</sup>–12</sup></td><td>[7.0, 12.0]</td><td>”</td><td>”</td><td>”</td><td>5–10 km</td></tr><tr><td>Config<sub>2</sub><sup>F<sup>7</sup>–12</sup>(3-D)</td><td>[7.0,12.0]</td><td>2 Myr</td><td>0.5–2 cm yr<sup>-1</sup></td><td>[700, 900; 100] °C</td><td>5–10 km</td></tr></tbody></table>

at the upper and lower edges of the convective domain, thus reducing gravitational contrasts and killing the buoyant instabilities faster.

Overall, we observe that in this Configuration 1, it is not sufficient to cool down the system to maintain the accumulated light inclusions ‘domes’ below the upper crust. We find it necessary to account for the ‘freezing’ of the inclusions density once cooling begins. This strong approximation is further discussed in Section 5.2.

### 4.1.4 Configuration 1 – Selected 3-D case

The selected configuration in 3-D is named *Config*₁<sup>8–13</sup>(3-D) (video via the link in Table 2). Initial and boundary conditions are similar to the 2-D cases. Here, basal heating decreases from $t_1 = 8$ Myr to $t_2 = 13$ Myr, and the light inclusions density freezes from time $t_1$.

Snapshots in time are displayed Fig. 5. The upper part of the figure presents the global geometry of the system in 3-D and its lower part displays 2-D profiles. At $t = 5.4$ Myr, convection has started and 3-D cells form, of horizontal size of *ca.* 20 km. These convective cells are responsible for the formation of large scale dome structures well visible at time $t = 8.1$ Myr. Then basal heating reduces, and the large scale domes keep evolving while the light inclusions start stacking at 15–20 km depth (at $t = 10.4$ and 14 Myr). In the later stages while the system cools down, the large scale domes vanish and the stacked small-scale domes maintain a rather heterogeneous texture with mixed light material

and BED material. The typical size of these structures is about 5 km.

Fig. 6 displays tracers trajectories in 3-D and their properties, showing the following features:

(i) The convective cycles periods are $\tau \approx 2$ Myr, lasting from about 5–12 Myr.

(ii) The local velocity magnitude is $2–3$ cm yr$^{-1}$ with a maximum of $10$ cm yr$^{-1}$ during the convecting period, reducing to a few mm yr$^{-1}$ after 13 Myr.

(iii) The tracers temperature difference over two cycles ranges between *ca.* 700 and 900 °C.

## 4.2 Configuration 2 – layering regime

We consider now Configuration 2, where the viscosity of the unmolten layer is chosen 2.3 times larger than in Configuration 1 (cf. parameters in Table 1), so that the system’s dynamics shifts towards the so-called layering regime identified in Louis-Napoléon et al. (2022). With such a smaller $Ra_{UM}$, we expect reduced convection vigour and velocities with respect to Configuration 1.

### 4.2.1 Configuration 2 – No cooling

This model case (*Config*₂<sup>∞</sup>) accounts for basal heating over the entire model duration, for example 20 Myr. Snapshots in time are displayed in Fig. 7(a). In contrast to Configuration 1, the light and heavy inclusions are able to escape from the convective cells and
<table><caption>**Table 1.** Parameters for Configurations 1 and 2: *i* is phase index (*i* = 1 for the BED, *i* = 2 for the light and less viscous inclusions and *i* = 3 for the heavy more viscous inclusions); PM, partially molten; UM, unmolten material; $\rho_i^{L0}$ and $\rho_i^{S0}$ are the reference densities used in eq. (10); $H_i^r$ is radioactive heat production, $\mu_i^0$ and $A_i$ are viscosity pre-factors in eqs (9) and (6), respectively (for Configurations 1 or 2, labelled in exponent). Furthermore, $n = 2.3$, $Q = 1.54 \times 10^5$ J mol$^{-1}$, $R = 8.314$ J mol$^{-1}$ K$^{-1}$, $\alpha = 3 \times 10^{-5}$ K$^{-1}$, so that $Ra_{UM} = 100$ and $Ra_{PM} = 3100$ in Configuration 1 and $Ra_{UM} = 13$ and $Ra_{PM} = 3100$ in Configuration 2.</caption><thead><tr><th>Fluid</th><th>Density<br>(PM layer)</th><th>Density<br>(UM layer)</th><th>Production</th><th>Viscosity<br>(PM layer)</th><th>Viscosity pre-factor<br>(UM layer)</th><th>Viscosity pre-factor<br>(UM layer)</th></tr><tr><th><i>i</i></th><th><i>ρ</i><sub><i>i</i></sub><sup>L0</sup><br>(kg m<sup>-3</sup>)</th><th><i>ρ</i><sub><i>i</i></sub><sup>S0</sup><br>(kg m<sup>-3</sup>)</th><th><i>H</i><sub><i>i</i></sub><sup>r</sup><br>(μW m<sup>-3</sup>)</th><th><i>μ</i><sub><i>i</i></sub><sup>0</sup><br>(Pa s)</th><th>Config<sub>1</sub>: <i>A</i><sub><i>i</i></sub><sup>1</sup><br>(Pa<sup>-2.3</sup> s<sup>-1</sup>)</th><th>Config<sub>2</sub>: <i>A</i><sub><i>i</i></sub><sup>2</sup><br>(Pa<sup>-2.3</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1</td><td>2700</td><td>2800</td><td>1.75</td><td>8.7 × 10<sup>16</sup></td><td>3.2 × 10<sup>-4</sup></td><td>4.2 × 10<sup>-5</sup></td></tr><tr><td>2</td><td>2500</td><td>2800</td><td>10</td><td>8.7 × 10<sup>15</sup></td><td>3.2 × 10<sup>-4</sup></td><td>4.2 × 10<sup>-5</sup></td></tr><tr><td>3</td><td>2900</td><td>2800</td><td>0.9</td><td>8.7 × 10<sup>17</sup></td><td>3.2 × 10<sup>-4</sup></td><td>4.2 × 10<sup>-5</sup></td></tr></tbody></table>

<table><caption>**Table 2.** Configurations names (*Config*₁ and *Config*₂) distinguish options: the exponents indicate the onset and ending times of cooling at the bottom boundary ($\infty$ = no cooling); ‘F’ refers to freezing of the inclusions density from cooling time $t_1$; (3-D) indicates a 3-D simulation. The models produce convection characterized by a ‘Cycle time’, for example the convection period, and ‘Velocity’, for example the mean markers velocity range. The two last columns display convective cells temperature range (prior to the onset of cooling and measured from a non-exhaustive selection of markers), and the final size of the domes when the system has come to rest.</caption><thead><tr><th></th><th>Simulation name</th><th>Cooling period<br>[t<sub>1</sub>, t<sub>2</sub>] (in Myr)</th><th>Cycle time (τ)</th><th>Convection properties<br>Velocity (v<sub>rms</sub>)</th><th>[T<sub>min</sub>, T<sub>max</sub>; ΔT]</th><th>Domes<br>D<sub>20Myr</sub></th></tr></thead><tbody><tr><td rowspan="5">Configuration 1</td><td>Config<sub>1</sub><sup>∞</sup></td><td>[∞, ∞]</td><td>1.8 Myr</td><td>2–5 cm yr<sup>-1</sup></td><td>[700, 850; 150] °C</td><td>Ø;</td></tr><tr><td>Config<sub>1</sub><sup>5–10</sup></td><td>[4.4,9.4]</td><td>”</td><td>”</td><td>”</td><td>Ø;</td></tr><tr><td>Config<sub>1</sub><sup>5–15</sup></td><td>[4.4,14.4]</td><td>”</td><td>”</td><td>”</td><td>Ø;</td></tr><tr><td>Config<sub>1</sub><sup>F<sup>5</sup>–10</sup></td><td>[4.4,9.4]</td><td>”</td><td>”</td><td>”</td><td>10 km</td></tr><tr><td>Config<sub>1</sub><sup>F<sup>8</sup>–13</sup>(3-D)</td><td>[8.0,13.0]</td><td>2 Myr</td><td>2–3 cm yr<sup>-1</sup></td><td>[700, 900; 200] °C</td><td>10 km</td></tr><tr><td rowspan="5">Configuration 2</td><td>Config<sub>2</sub><sup>∞</sup></td><td>[∞, ∞]</td><td>2 Myr</td><td>0.5–1 cm yr<sup>-1</sup></td><td>[700, 850; 150] °C</td><td>5–10 km</td></tr><tr><td>Config<sub>2</sub><sup>7–12</sup></td><td>[7.0, 12.0]</td><td>”</td><td>”</td><td>”</td><td>5–10 km</td></tr><tr><td>Config<sub>2</sub><sup>F<sup>7</sup>–12</sup></td><td>[7.0, 12.0]</td><td>”</td><td>”</td><td>”</td><td>5–10 km</td></tr><tr><td>Config<sub>2</sub><sup>F<sup>7</sup>–12</sup>(3-D)</td><td>[7.0,12.0]</td><td>2 Myr</td><td>0.5–2 cm yr<sup>-1</sup></td><td>[700, 900; 100] °C</td><td>5–10 km</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

at the upper and lower edges of the convective domain, thus reducing gravitational contrasts and killing the buoyant instabilities faster.

Overall, we observe that in this Configuration 1, it is not sufficient to cool down the system to maintain the accumulated light inclusions ‘domes’ below the upper crust. We find it necessary to account for the ‘freezing’ of the inclusions density once cooling begins. This strong approximation is further discussed in Section 5.2.

### 4.1.4 Configuration 1 – Selected 3-D case

The selected configuration in 3-D is named *Config*₁<sup>8–13</sup>(3-D) (video via the link in Table 2). Initial and boundary conditions are similar to the 2-D cases. Here, basal heating decreases from $t_1 = 8$ Myr to $t_2 = 13$ Myr, and the light inclusions density freezes from time $t_1$.

Snapshots in time are displayed Fig. 5. The upper part of the figure presents the global geometry of the system in 3-D and its lower part displays 2-D profiles. At $t = 5.4$ Myr, convection has started and 3-D cells form, of horizontal size of *ca.* 20 km. These convective cells are responsible for the formation of large scale dome structures well visible at time $t = 8.1$ Myr. Then basal heating reduces, and the large scale domes keep evolving while the light inclusions start stacking at 15–20 km depth (at $t = 10.4$ and 14 Myr). In the later stages while the system cools down, the large scale domes vanish and the stacked small-scale domes maintain a rather heterogeneous texture with mixed light material

and BED material. The typical size of these structures is about 5 km.

Fig. 6 displays tracers trajectories in 3-D and their properties, showing the following features:

(i) The convective cycles periods are $\tau \approx 2$ Myr, lasting from about 5–12 Myr.

(ii) The local velocity magnitude is $2–3$ cm yr$^{-1}$ with a maximum of $10$ cm yr$^{-1}$ during the convecting period, reducing to a few mm yr$^{-1}$ after 13 Myr.

(iii) The tracers temperature difference over two cycles ranges between *ca.* 700 and 900 °C.

## 4.2 Configuration 2 – layering regime

We consider now Configuration 2, where the viscosity of the unmolten layer is chosen 2.3 times larger than in Configuration 1 (cf. parameters in Table 1), so that the system’s dynamics shifts towards the so-called layering regime identified in Louis-Napoléon et al. (2022). With such a smaller $Ra_{UM}$, we expect reduced convection vigour and velocities with respect to Configuration 1.

### 4.2.1 Configuration 2 – No cooling

This model case (*Config*₂<sup>∞</sup>) accounts for basal heating over the entire model duration, for example 20 Myr. Snapshots in time are displayed in Fig. 7(a). In contrast to Configuration 1, the light and heavy inclusions are able to escape from the convective cells and
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Configuration 1, 3-D model *Config1F*⁸⁻¹³ (3-D). Top rows: six snapshots in time in 3-D view display light less viscous (white) and heavy more viscous (black) inclusions, and isotherm *T* = 720 °C (iso-surface associated to the melting front). Bottom rows: snapshots in 2-D set in the *z* = 25 km plane, with temperature isocontours, the viscosity (left half of model domain), and the light (white) and heavy (black) inclusions (right half of the model domain).</figcaption></figure>

form layers below the upper crust and at the bottom of the domain, respectively. This behaviour is typical of the layering regime (Louis-Napoléon *et al.* 2022). A closer look at the temporal evolution of the layer formed by the light inclusions reveals that it is strongly deformed by the convective flow. Over time, this layer adopts a flat shape, destroying the dome structures.

### 4.2.2 Configuration 2 – With cooling

Now, we progressively reduce the basal heating, starting from time *t*₁ = 7 Myr until time *t*₂ = 12 Myr (Fig. 7b). This cooling time is delayed with respect to Configuration 1 (where *t*₁ = 4.4 Myr) because we need to let time for a few convection cycles to develop,
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Configuration 1, 3-D model *Config1F*⁸⁻¹³ (3-D). Top rows: six snapshots in time in 3-D view display light less viscous (white) and heavy more viscous (black) inclusions, and isotherm *T* = 720 °C (iso-surface associated to the melting front). Bottom rows: snapshots in 2-D set in the *z* = 25 km plane, with temperature isocontours, the viscosity (left half of model domain), and the light (white) and heavy (black) inclusions (right half of the model domain).</figcaption></figure>

form layers below the upper crust and at the bottom of the domain, respectively. This behaviour is typical of the layering regime (Louis-Napoléon *et al.* 2022). A closer look at the temporal evolution of the layer formed by the light inclusions reveals that it is strongly deformed by the convective flow. Over time, this layer adopts a flat shape, destroying the dome structures.

### 4.2.2 Configuration 2 – With cooling

Now, we progressively reduce the basal heating, starting from time *t*₁ = 7 Myr until time *t*₂ = 12 Myr (Fig. 7b). This cooling time is delayed with respect to Configuration 1 (where *t*₁ = 4.4 Myr) because we need to let time for a few convection cycles to develop,
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Configuration 1: Time evolution of some markers in the 3-D model *Config1F*<sup>8</sup>–<sup>13</sup>(3-D). Their initial location is marked by circles in frame (a). (a) 3-D trajectory, (b) depth, (c) temperature and (d) local velocity magnitude. Here, the convection revolving period is about 2 Myr, with mean velocities of 2–3 cm yr<sup>-1</sup>, and maximum velocities of 10 cm yr<sup>-1</sup>. Temperatures vary by almost 200 °C around about 800 °C.</figcaption></figure>

else convection is inhibited. As expected, thermal convection pro-
gressively dies out and the layer of light inclusions formed at 18 km
depth becomes less deformed over time. Hence, irregular structures
of light inclusions of characteristic sizes *ca*. 5 km are preserved, but
form neither flat layers nor dome structures.

4.2.3 Configuration 2 – With cooling and freezing properties of inclusions

We impose now that the density of the white inclusions is maintained
fixed from t₁, as tested in Configuration 1. Results are displayed in
Fig. 7(c). From time ca. 15 Myr, domes become clearly visible,
and reach about 18 km depth at *t* = 19.3 Myr. Note also that
some domains of the BED fluid remain trapped within the domes,
illustrating the preservation of the system’s heterogeneity. Here
again, the typical size of the domes is 5 km.

4.2.4 Configuration 2 – Selected 3-D case

The selected configuration in 3-D is named *Config2F*<sup>7–12</sup>(3-D) (video via the link in Table 2). Initial and boundary conditions are similar to the 2-D case, with basal cooling from *t*₁ ≈ 7 Myr until *t*₂ ≈ 12 Myr, and freezing of the inclusions density from time *t*₁.

Snapshots in time are displayed in Fig. 8, with 3-D views in its
upper part and 2-D vertical profiles in its lower part. In this layering
regime, the convection vigour is significantly reduced compared to
that in the suspension regime (Configuration 1). This can be seen
when comparing the mean vertical velocity between both configu-
rations (see e.g. Figs 6d and 9d), which differ by almost a factor
2, as well as when comparing the viscosity ranges, for instance
around 20 km depth. Here, because convection initiates later, the
light inclusions segregate and cluster under the melting front prior
to 5 Myr, and these clusters display a somewhat rounded shape of
about 2–5 km in size (cf. 2-D profiles at 5 and 6 Myr), despite the
development of convection. After 14 Myr and as the system cools
down, heterogeneous structures of sizes 2–5 km still lie at 18–22 km
depth.

Fig. 9 displays some markers trajectories in 3-D view and their
properties, which allow to deduce the following information:

(i) The period of the convective cycles appears irregular ranging
in $\tau \approx 1-2$ Myr, and lasts from about 5–12 Myr.

(ii) The local velocity magnitude is 0.5–2 cm yr⁻¹ with a maximum at 5 cm yr⁻¹, reducing to a few mm yr⁻¹ after 12 Myr.

(iii) The tracers temperature during three cycles varies between 700 and 900 °C, but the amplitude from one cycle to the other reaches about 100 °C.

5 SYNTHESIS AND DISCUSSION

**5.1 Models synthesis**

The models presented above demonstrate the possibility for par-
tially molten crust material to stack in the form of domes at depths
around 15 km, provided an appropriate range of rheological and
boundary conditions. In general, these domes have a size of approx-
imately 5–10 km, they form and stabilize within about 10 Myr, and
they record several thermal peaks which reflect convective motion.
Such convecting cycles provide an explanation for the cyclic pre-
cipitation/dissolution record of sampled zircons on Naxos Island
(Vanderhaeghe *et al.* 2018).

Accounting for crustal convection is necessary if one wishes
to explain the observation of low thermal gradients over broad
migmatite domains, as mentioned in Section 2. For instance let us
recall Riel *et al.* (2016) conclusions for the El Oro Triassic metased-
imentary complex in Ecuador: pressure and temperature estimates
show a gradient of *ca*. 10 °C km⁻¹ in the 7–8 km thick garnet-
bearing migmatitic lower crust. According to petrological and geo-
chemical indicators, melts produced during biotite-breakdown (5–
15 vol. per cent melt) were trapped and pervasively distributed in
these rocks. Riel *et al.* (2016) could not explain such a low thermal
gradient with neither diffusive nor upward melt transfer processes
and hence concluded that crustal convection must have occurred.
Furthermore, Vanderhaeghe *et al.* (2018) provided analytical esti-
mates of the minimum viscosity required for convection to develop
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Configuration 1: Time evolution of some markers in the 3-D model *Config1F*<sup>8</sup>–<sup>13</sup>(3-D). Their initial location is marked by circles in frame (a). (a) 3-D trajectory, (b) depth, (c) temperature and (d) local velocity magnitude. Here, the convection revolving period is about 2 Myr, with mean velocities of 2–3 cm yr<sup>-1</sup>, and maximum velocities of 10 cm yr<sup>-1</sup>. Temperatures vary by almost 200 °C around about 800 °C.</figcaption></figure>

else convection is inhibited. As expected, thermal convection pro-
gressively dies out and the layer of light inclusions formed at 18 km
depth becomes less deformed over time. Hence, irregular structures
of light inclusions of characteristic sizes *ca*. 5 km are preserved, but
form neither flat layers nor dome structures.

4.2.3 Configuration 2 – With cooling and freezing properties of inclusions

We impose now that the density of the white inclusions is maintained
fixed from t₁, as tested in Configuration 1. Results are displayed in
Fig. 7(c). From time ca. 15 Myr, domes become clearly visible,
and reach about 18 km depth at *t* = 19.3 Myr. Note also that
some domains of the BED fluid remain trapped within the domes,
illustrating the preservation of the system’s heterogeneity. Here
again, the typical size of the domes is 5 km.

4.2.4 Configuration 2 – Selected 3-D case

The selected configuration in 3-D is named *Config2F*<sup>7–12</sup>(3-D) (video via the link in Table 2). Initial and boundary conditions are similar to the 2-D case, with basal cooling from *t*₁ ≈ 7 Myr until *t*₂ ≈ 12 Myr, and freezing of the inclusions density from time *t*₁.

Snapshots in time are displayed in Fig. 8, with 3-D views in its
upper part and 2-D vertical profiles in its lower part. In this layering
regime, the convection vigour is significantly reduced compared to
that in the suspension regime (Configuration 1). This can be seen
when comparing the mean vertical velocity between both configu-
rations (see e.g. Figs 6d and 9d), which differ by almost a factor
2, as well as when comparing the viscosity ranges, for instance
around 20 km depth. Here, because convection initiates later, the
light inclusions segregate and cluster under the melting front prior
to 5 Myr, and these clusters display a somewhat rounded shape of
about 2–5 km in size (cf. 2-D profiles at 5 and 6 Myr), despite the
development of convection. After 14 Myr and as the system cools
down, heterogeneous structures of sizes 2–5 km still lie at 18–22 km
depth.

Fig. 9 displays some markers trajectories in 3-D view and their
properties, which allow to deduce the following information:

(i) The period of the convective cycles appears irregular ranging
in $\tau \approx 1-2$ Myr, and lasts from about 5–12 Myr.

(ii) The local velocity magnitude is 0.5–2 cm yr⁻¹ with a maximum at 5 cm yr⁻¹, reducing to a few mm yr⁻¹ after 12 Myr.

(iii) The tracers temperature during three cycles varies between 700 and 900 °C, but the amplitude from one cycle to the other reaches about 100 °C.

5 SYNTHESIS AND DISCUSSION

**5.1 Models synthesis**

The models presented above demonstrate the possibility for par-
tially molten crust material to stack in the form of domes at depths
around 15 km, provided an appropriate range of rheological and
boundary conditions. In general, these domes have a size of approx-
imately 5–10 km, they form and stabilize within about 10 Myr, and
they record several thermal peaks which reflect convective motion.
Such convecting cycles provide an explanation for the cyclic pre-
cipitation/dissolution record of sampled zircons on Naxos Island
(Vanderhaeghe et al. 2018).

Accounting for crustal convection is necessary if one wishes
to explain the observation of low thermal gradients over broad
migmatite domains, as mentioned in Section 2. For instance let us
recall Riel *et al.* (2016) conclusions for the El Oro Triassic metased-
imentary complex in Ecuador: pressure and temperature estimates
show a gradient of *ca*. 10 °C km⁻¹ in the 7–8 km thick garnet-
bearing migmatitic lower crust. According to petrological and geo-
chemical indicators, melts produced during biotite-breakdown (5–
15 vol. per cent melt) were trapped and pervasively distributed in
these rocks. Riel *et al.* (2016) could not explain such a low thermal
gradient with neither diffusive nor upward melt transfer processes
and hence concluded that crustal convection must have occurred.
Furthermore, Vanderhaeghe *et al.* (2018) provided analytical esti-
mates of the minimum viscosity required for convection to develop
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Configuration 2: Summary of the 2-D results. From top to bottom: (a) *Config<sub>2</sub>*<sup>∞</sup>, (b) *Config<sub>2</sub>*<sup>7−12</sup> and (c) *Config<sub>2</sub>*F<sup>7−12</sup> (see Tables 1 and 2 for parameters). For each case, snapshots display temperature isocontours, in the left half of the model domain the viscosity distribution, and, in the right half, the light and heavy inclusions in white and black, respectively. (d) Time evolution of the average velocity.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Configuration 2: Summary of the 2-D results. From top to bottom: (a) *Config<sub>2</sub>*<sup>∞</sup>, (b) *Config<sub>2</sub>*<sup>7−12</sup> and (c) *Config<sub>2</sub>*F<sup>7−12</sup> (see Tables 1 and 2 for parameters). For each case, snapshots display temperature isocontours, in the left half of the model domain the viscosity distribution, and, in the right half, the light and heavy inclusions in white and black, respectively. (d) Time evolution of the average velocity.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Configuration 2: Time evolution of the flow in the 3-D model *Config2F*<sup>7</sup>–<sup>12</sup>(3-D) ('layering' regime, basal heating and cooling, and freezing inclusions density), in 3-D (top panels) and 2-D vertical sections (bottom panels). Same legend as Fig. 5.</figcaption></figure>

over a range of thicknesses of the Naxos crust (Greece). Our models here show that convection can indeed develop over thicknesses up to about 20 km in the lower crust within about 10 Myr of basal heating.

The two configurations tested above also display differences in terms of domes dynamics. We summarize here the main differences,

considering the 3-D models (which account for basal cooling and freezing of inclusions density):

(i) In the suspension regime (Configuration 1), the domes have sizes of about 10 km. Cyclic revolution records thermal oscillations around 700–900 °C over periods of *ca.* 1.5 Myr, starting 4–5 Myr after the onset of basal heating. Vertical rms velocities are of the
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Configuration 2: Time evolution of the flow in the 3-D model *Config2F*<sup>7-12</sup>(3-D) ('layering' regime, basal heating and cooling, and freezing inclusions density), in 3-D (top panels) and 2-D vertical sections (bottom panels). Same legend as Fig. 5.</figcaption></figure>

over a range of thicknesses of the Naxos crust (Greece). Our models here show that convection can indeed develop over thicknesses up to about 20 km in the lower crust within about 10 Myr of basal heating.

The two configurations tested above also display differences in terms of domes dynamics. We summarize here the main differences,

considering the 3-D models (which account for basal cooling and freezing of inclusions density):

(i) In the suspension regime (Configuration 1), the domes have sizes of about 10 km. Cyclic revolution records thermal oscillations around 700–900 °C over periods of *ca.* 1.5 Myr, starting 4–5 Myr after the onset of basal heating. Vertical rms velocities are of the
<figure><img src="image_9.png" /><figcaption>**Figure 9.** Configuration 2: Time evolution of three markers in the selected 3-D model *Config2F*<sup>7</sup>–<sup>12</sup>(3-D). Their initial location is marked by circles in frame (a). (a) 3-D trajectory, (b) depth location, (c) temperature and (d) local velocity magnitude. The convective revolving period is about 2 Myr, at a velocity in the range 0.5–2 cm·yr<sup>-1</sup>, maximum velocities reach 5 cm·yr<sup>-1</sup>, and temperatures vary by 100 °C in the range [700, 900] °C. Convection ceases at time *t* ≈ 12 Myr.</figcaption></figure>

order of 2–4 cm yr⁻¹. The domes shape is formed by aggregated light inclusions with enclosed pieces of the BED with layers of different properties. They adopt their final shape from time *ca.* 10 Myr and reach an exhumation depth of *ca.* 12–15 km after 20 Myr.

(ii) In the layering regime (Configuration 2), domes have smaller sizes of about 2–5 km, similar thermal variations are slightly lower but occur over a similar range around 800 °C, starting at ~5 Myr after the onset of basal heating. Vertical rms velocities stand in the range of 0.5–2 cm yr⁻¹. The domes have a more homogeneous composition, and they reach depth 15–18 km after 20 Myr.

The suspension regime is likely to correspond to orogenic contexts in which the crust heats up enough to undergo vigorous thermal convection. Such a configuration can be associated with a generally low viscosity behaviour and rather high Rayleigh numbers, that applies rather well to ancient geodynamic contexts. The formation of dome structures is favoured when accounting for the following additional processes: (i) an appropriate timing of the subsequent basal cooling and (ii) sufficient melt extraction at the inclusions subscale so that the inclusions density becomes ‘melt-insensitive’ at the several Myr time-scale. Hence the modelled domes, represented by a high concentration of white inclusions, should be seen as a heterogeneous mixture between initial crustal compositions and partially molten material.

The layering regime produces smaller dome structures, and is characterized by generally greater viscosities and lower Rayleigh numbers than Configuration 1. The dynamics are significantly slower (e.g. mm yr⁻¹ as opposed to cm yr⁻¹ in the suspension regime), which better corresponds to present-day geodynamic contexts; such modern partially molten crustal roots would still have time to mature thermally prior to an external tectonic forcing (e.g. lateral extension). These domes can remain stuck at mid-crustal level for some time before lateral extensional tectonics drive further exhumation to the surface.

Finally, 3-D models are necessary to consider since they capture the inherent 3-D dynamics and rounded structure of the compositional domes, as opposed to the 2-D models which can only reproduce 2-D tubular-like structures. Configuration 1 reproduces the

formation of nested small domes into a larger upwelling dome (a convection cell). Furthermore, in both configurations, convection and dome structures display 3-D ‘radial’ patterns that are comparable to the stretching lineations observed in natural examples (cf. citations in Section 2 and cf. Louis-Napoléon 2020).

## 5.2 Limits of the modelling approximations

Our models make a number of assumptions which deserve discussion, prior to extrapolation to real geodynamic systems. Some parameters choices have not been discussed in detail because first, a large number of individual parameters have already been varied and tackled in previous studies (see below); secondly, for many of these parameters a precise value is not achievable given the variety and non-linearity of felsic rocks behaviour; third, we are looking for key trends rather than detailed values that would not bring further insight on the major controls on migmatite dome formation. Nonetheless, we have tested the influence of a number of parameters : for instance the conditions for polydiapirism and the influence of a third dimension on Rayleigh–Taylor and Rayleigh–Bénard instabilities were tested in Louis-Napoléon et al. (2020). The distribution of the inclusions within the model domain and a variety of viscosities were tested in Louis-Napoléon et al. (2022). Let us address below other important factors dealing with the inclusions properties: melt fraction dependent rheology and the inclusions concentration. We also discuss the role of the boundary conditions.

### 5.2.1 Other parameters that determine the suspension and layering regimes

In this study, we selected two configurations that only differ by their ‘unmolten’ viscosity so that one falls in the suspension regime and the other in the layering regime. Configuration 1 has smaller unmolten viscosity hence thinner stiff lid than Configuration 2; consequently, convection occurs over a thicker layer, the effective Ra<sub>UM</sub> is higher and the rms velocity higher than in Configuration 2. This leads to relatively slower segregation of light inclusions beneath the upper crustal lid, with respect to Configuration 2. In
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Configuration 2: Time evolution of three markers in the selected 3-D model *Config2F*<sup>7</sup>–<sup>12</sup>(3-D). Their initial location is marked by circles in frame (a). (a) 3-D trajectory, (b) depth location, (c) temperature and (d) local velocity magnitude. The convective revolving period is about 2 Myr, at a velocity in the range 0.5–2 cm·yr<sup>−1</sup>, maximum velocities reach 5 cm·yr<sup>−1</sup>, and temperatures vary by 100 °C in the range [700, 900] °C. Convection ceases at time *t* ≈ 12 Myr.</figcaption></figure>

order of 2–4 cm yr⁻¹. The domes shape is formed by aggregated light inclusions with enclosed pieces of the BED with layers of different properties. They adopt their final shape from time *ca*. 10 Myr and reach an exhumation depth of *ca*. 12–15 km after 20 Myr.

(ii) In the layering regime (Configuration 2), domes have smaller sizes of about 2–5 km, similar thermal variations are slightly lower but occur over a similar range around 800 °C, starting at ~5 Myr after the onset of basal heating. Vertical rms velocities stand in the range of 0.5–2 cm yr⁻¹. The domes have a more homogeneous composition, and they reach depth 15–18 km after 20 Myr.

The suspension regime is likely to correspond to orogenic contexts in which the crust heats up enough to undergo vigorous thermal convection. Such a configuration can be associated with a generally low viscosity behaviour and rather high Rayleigh numbers, that applies rather well to ancient geodynamic contexts. The formation of dome structures is favoured when accounting for the following additional processes: (i) an appropriate timing of the subsequent basal cooling and (ii) sufficient melt extraction at the inclusions subscale so that the inclusions density becomes ‘melt-insensitive’ at the several Myr time-scale. Hence the modelled domes, represented by a high concentration of white inclusions, should be seen as a heterogeneous mixture between initial crustal compositions and partially molten material.

The layering regime produces smaller dome structures, and is characterized by generally greater viscosities and lower Rayleigh numbers than Configuration 1. The dynamics are significantly slower (e.g. mm yr⁻¹ as opposed to cm yr⁻¹ in the suspension regime), which better corresponds to present-day geodynamic contexts; such modern partially molten crustal roots would still have time to mature thermally prior to an external tectonic forcing (e.g. lateral extension). These domes can remain stuck at mid-crustal level for some time before lateral extensional tectonics drive further exhumation to the surface.

Finally, 3-D models are necessary to consider since they capture the inherent 3-D dynamics and rounded structure of the compositional domes, as opposed to the 2-D models which can only reproduce 2-D tubular-like structures. Configuration 1 reproduces the

formation of nested small domes into a larger upwelling dome (a convection cell). Furthermore, in both configurations, convection and dome structures display 3-D ‘radial’ patterns that are comparable to the stretching lineations observed in natural examples (cf. citations in Section 2 and cf. Louis-Napoléon 2020).

## 5.2 Limits of the modelling approximations

Our models make a number of assumptions which deserve discussion, prior to extrapolation to real geodynamic systems. Some parameters choices have not been discussed in detail because first, a large number of individual parameters have already been varied and tackled in previous studies (see below); secondly, for many of these parameters a precise value is not achievable given the variety and non-linearity of felsic rocks behaviour; third, we are looking for key trends rather than detailed values that would not bring further insight on the major controls on migmatite dome formation. Nonetheless, we have tested the influence of a number of parameters : for instance the conditions for polydiapirism and the influence of a third dimension on Rayleigh–Taylor and Rayleigh–Bénard instabilities were tested in Louis-Napoléon et al. (2020). The distribution of the inclusions within the model domain and a variety of viscosities were tested in Louis-Napoléon et al. (2022). Let us address below other important factors dealing with the inclusions properties: melt fraction dependent rheology and the inclusions concentration. We also discuss the role of the boundary conditions.

### 5.2.1 Other parameters that determine the suspension and layering regimes

In this study, we selected two configurations that only differ by their ‘unmolten’ viscosity so that one falls in the suspension regime and the other in the layering regime. Configuration 1 has smaller unmolten viscosity hence thinner stiff lid than Configuration 2; consequently, convection occurs over a thicker layer, the effective Ra<sub>UM</sub> is higher and the rms velocity higher than in Configuration 2. This leads to relatively slower segregation of light inclusions beneath the upper crustal lid, with respect to Configuration 2. In
paper part I (Louis-Napoléon et al. 2022), we explored a number of other parameters that led us to the regime diagram of Fig. 3. Let us recall below the other key parameters that control the system’s evolution :

(i) The Stokes number $St$, defined as the ratio of the characteristic response time of the inclusions $\tau_{incl}$ to the characteristic time of the convective flow $\tau_{flow}$, influences the system’s evolution. If $St \ll 1$, inclusions are likely to follow the local flow, whereas if $St \gg 1$ they will move independently. As a rough estimate, we can write $St = \tau_{incl}/\tau_{flow}$ with $\tau_{incl} = 2\rho r^2/(9\mu)$ and $\tau_{flow} = (H_T/2)/v_{flow}$, where $v_{flow} = \sqrt{\alpha g \Delta T (H_T/2)}$ characterizes the convective flow velocity. Using $\rho = 2700 \text{ kg m}^{-3}$, $r = 600 \text{ m}$, $\mu = 10^{17} \text{ Pa s}$, $H_T = 45 \text{ km}$, $\alpha = 3 \times 10^{-5} \text{ K}^{-1}$, $g = 9.81 \text{ m s}^{-2}$, $\Delta T = 400 \text{ K}$, we obtain $St \approx 5 \times 10^{-12}$. This shows that the inclusions inertia is negligible here and cannot be responsible for the segregation or layering processes. Moreover, we can compute the inclusions characteristic Stokes velocity $v_s = \frac{2}{9} \frac{\Delta \rho g r^2}{\mu}$ and characteristic residence time $\tau_s = (H_T/2)/v_s$. With $\Delta \rho = 200 \text{ kg m}^{-3}$ and the values given above, we obtain $v_s \approx 5 \text{ cm yr}^{-1}$ and $\tau_s \approx 0.5 \text{ Myr}$. We showed in Louis-Napoléon et al. (2022) that in our setting of convective flow, with $St \approx 5 \times 10^{-12}$ and a buoyancy ratio $\Lambda < 10$, our models stand in the so called ‘dust-like’ regime according to Patočka et al. (2020); in this regime the inclusions tend to segregate within a characteristic time $3\tau_s$, which approaches 1.5 Myr here. Since this time duration is shorter than the time over which convection develops in our models (about 4 Myr), the inclusions have time to segregate and/or form layers.

(ii) The viscosity and density ratios between the inclusions and the surrounding material (BED) control whether the inclusions are able to segregate into layers efficiently. The sensitivity of these parameters was only superficially explored in Part I. Investigating the influence of viscosity and density ratios in more detail is needed, but stands beyond the scope of this study.

(iii) The intensity of convection is affected by the volume of buoyant melt generated. This depends on the solidus and liquidus temperatures. $T_{sol}$ and $T_{liq}$ are linearly linked to $M$ and $M_s$ and intervene in the Rayleigh numbers via the melting front temperature $T(M = M_s)$. Hence varying $T_{sol}$ and $T_{liq}$ is equivalent to varying $M_s$ and $R_a$. In fact Appendix A1.1 in part I (Louis-Napoléon et al. 2022), displays a sensibility analysis on $M_s$ for $M_s = 0.1$ and 0.3 (see case $I_0 M_1 H_0$ versus $I_0 M_1 H_5$ in Table A2), and shows that the main behaviours remain. However it would be useful to test phase dependent melting rates ($M(T_{liq}^i, T_{sol}^i)$).

All these parameters influence the boundaries between the system’s flow regimes (cf. Fig. 3 and convection onset time. Although here we have only investigated the conditions for dome formation by simply varying the viscosity as a switch between the suspension and layering regimes, it is clear that other parameters impacting the inclusions density, viscosity, as well as the thermal state, also modulate this transition from one regime to the other.

## 5.2.2 Melt properties – influence of the inclusions size and concentration

We have explained our strong assumption of how our modelled inclusions are supposed to capture melt percolation processes at the scale greater than several hundred metres. Although we relied on the sizes predicted by previous studies (cf. Section 3), further work requires to better constrain the conditions by which this assumption holds. Although the threshold $M_s = 0.1$ or 0.3 was shown not to

drastically change the system’s dynamics (Louis-Napoléon et al. 2022), this ‘structure’ of melt heterogeneity obviously requires an appropriate amount of felsic and hydrated minerals with respect to mafic components, to ensure not only sufficiently low overall viscosity, but also sufficient fluid pressure and density contrasts for melt segregation processes at the inclusions subscale (via porosity waves, e.g. Räss et al. 2019) or continuous percolation (Deubelbeiss & Kaus 2008; Brown 2010; Parmigiani et al. 2014). Several parameters actually have competing influences and addressing them here escapes the scope of this contribution. Nevertheless, future studies should improve our understanding of these intermediate scale processes.

We also tested how the concentration and size of the inclusions influences the layering regime and the formation of a layer of buoyant material at ca. 15 km depth, in Appendix C. Below a certain amount of inclusions, their presence does not influence the general flow and segregation dynamics, whereas above a maximum amount, they block the flow and impede segregation (consistent with e.g. Lavorel & Le Bars 2009; Harada et al. 2012). In Appendix C, we explain that an increase in the inclusions concentration $\phi$ tends to reduce the thickness of the convecting layer. Indeed, settling of the heavy inclusions occurs over a thicker layer (Fig. C1) while the melting front remains deeper (Fig. C2): both effects reduce the thickness of the convective layer. This in turn decreases the Rayleigh number of the unmolten region, reducing convection intensity and favouring the layering of light inclusions.

## 5.2.3 The crust's heterogeneity upon melting

Our Configurations 1 and 2 assume a relatively complex dependency of the rheology on melt fraction. However, such a switch in rheology from a temperature and strain-rate power-law relationship (eq. 6) to a melt fraction dependency (eq. 9) has become common practice in numerical geodynamics (e.g. Gerya et al. 2008; Perchuk et al. 2011; Schenker et al. 2012; Gerbault et al. 2018). Here we showed the key role of the crust’s heterogeneity in the formation of migmatite domes.

It would be interesting to carry a comparative study with other modelling approaches of partial melting continental crust. Such studies usually seek other purposes, such as mobile belt geodynamics in Archean and Proterozoic times (Sizova et al. 2015; Piccolo et al. 2019) or high pressure exhumation patterns in continental collision belts (Perchuk et al. 2011; Burov et al. 2014). Other theoretical models investigate multi-scale melt segregation processes (Räss et al. 2019; Schmeling et al. 2019; Culha et al. 2020; Schmeling et al. 2023) with advanced computational tools in order to deal with high resolution scales. It would be useful to investigate how these approaches also reproduce the segregation of partially molten felsic domes above a convective crust, in comparison to field estimates of dimensions and timing associated to migmatite terranes (e.g. Ledru et al. 2001; Riel et al. 2016; Vanderhaeghe et al. 2018). Such a comparison will help to confirm the requirements obtained in our study, regarding both the compositional heterogeneities and the boundary conditions.

## 5.2.4 The influence of the boundary conditions

The lateral boundaries were kept fixed in our models, consistent with the assumption that the orogenic system is in a transitional stage between continental convergence and post-orogenic collapse. Actually such a transitional stage may last several tens of millions
paper part I (Louis-Napoléon et al. 2022), we explored a number of other parameters that led us to the regime diagram of Fig. 3. Let us recall below the other key parameters that control the system’s evolution :

(i) The Stokes number $St$, defined as the ratio of the characteristic response time of the inclusions $\tau_{incl}$ to the characteristic time of the convective flow $\tau_{flow}$, influences the system’s evolution. If $St \ll 1$, inclusions are likely to follow the local flow, whereas if $St \gg 1$ they will move independently. As a rough estimate, we can write $St = \tau_{incl}/\tau_{flow}$ with $\tau_{incl} = 2\rho r^2/(9\mu)$ and $\tau_{flow} = (H_T/2)/v_{flow}$, where $v_{flow} = \sqrt{\alpha g \Delta T (H_T/2)}$ characterizes the convective flow velocity. Using $\rho = 2700 \text{ kg m}^{-3}$, $r = 600 \text{ m}$, $\mu = 10^{17} \text{ Pa s}$, $H_T = 45 \text{ km}$, $\alpha = 3 \times 10^{-5} \text{ K}^{-1}$, $g = 9.81 \text{ m s}^{-2}$, $\Delta T = 400 \text{ K}$, we obtain $St \approx 5 \times 10^{-12}$. This shows that the inclusions inertia is negligible here and cannot be responsible for the segregation or layering processes. Moreover, we can compute the inclusions characteristic Stokes velocity $v_s = \frac{2}{9} \frac{\Delta \rho g r^2}{\mu}$ and characteristic residence time $\tau_s = (H_T/2)/v_s$. With $\Delta \rho = 200 \text{ kg m}^{-3}$ and the values given above, we obtain $v_s \approx 5 \text{ cm yr}^{-1}$ and $\tau_s \approx 0.5 \text{ Myr}$. We showed in Louis-Napoléon et al. (2022) that in our setting of convective flow, with $St \approx 5 \times 10^{-12}$ and a buoyancy ratio $\Lambda < 10$, our models stand in the so called ‘dust-like’ regime according to Patočka et al. (2020); in this regime the inclusions tend to segregate within a characteristic time $3\tau_s$, which approaches 1.5 Myr here. Since this time duration is shorter than the time over which convection develops in our models (about 4 Myr), the inclusions have time to segregate and/or form layers.

(ii) The viscosity and density ratios between the inclusions and the surrounding material (BED) control whether the inclusions are able to segregate into layers efficiently. The sensitivity of these parameters was only superficially explored in Part I. Investigating the influence of viscosity and density ratios in more detail is needed, but stands beyond the scope of this study.

(iii) The intensity of convection is affected by the volume of buoyant melt generated. This depends on the solidus and liquidus temperatures. $T_{sol}$ and $T_{liq}$ are linearly linked to $M$ and $M_s$ and intervene in the Rayleigh numbers via the melting front temperature $T(M = M_s)$. Hence varying $T_{sol}$ and $T_{liq}$ is equivalent to varying $M_s$ and $R_a$. In fact Appendix A1.1 in part I (Louis-Napoléon et al. 2022), displays a sensibility analysis on $M_s$ for $M_s = 0.1$ and $0.3$ (see case $I_0 M_1 H_0$ versus $I_0 M_1 H_5$ in Table A2), and shows that the main behaviours remain. However it would be useful to test phase dependent melting rates ($M(T_{liq}^i, T_{sol}^i)$).

All these parameters influence the boundaries between the system’s flow regimes (cf. Fig. 3 and convection onset time. Although here we have only investigated the conditions for dome formation by simply varying the viscosity as a switch between the suspension and layering regimes, it is clear that other parameters impacting the inclusions density, viscosity, as well as the thermal state, also modulate this transition from one regime to the other.

## 5.2.2 Melt properties – influence of the inclusions size and concentration

We have explained our strong assumption of how our modelled inclusions are supposed to capture melt percolation processes at the scale greater than several hundred metres. Although we relied on the sizes predicted by previous studies (cf. Section 3), further work requires to better constrain the conditions by which this assumption holds. Although the threshold $M_s = 0.1$ or $0.3$ was shown not to

drastically change the system’s dynamics (Louis-Napoléon et al. 2022), this ‘structure’ of melt heterogeneity obviously requires an appropriate amount of felsic and hydrated minerals with respect to mafic components, to ensure not only sufficiently low overall viscosity, but also sufficient fluid pressure and density contrasts for melt segregation processes at the inclusions subscale (via porosity waves, e.g. Räss et al. 2019) or continuous percolation (Deubelbeiss & Kaus 2008; Brown 2010; Parmigiani et al. 2014). Several parameters actually have competing influences and addressing them here escapes the scope of this contribution. Nevertheless, future studies should improve our understanding of these intermediate scale processes.

We also tested how the concentration and size of the inclusions influences the layering regime and the formation of a layer of buoyant material at ca. 15 km depth, in Appendix C. Below a certain amount of inclusions, their presence does not influence the general flow and segregation dynamics, whereas above a maximum amount, they block the flow and impede segregation (consistent with e.g. Lavorel & Le Bars 2009; Harada et al. 2012). In Appendix C, we explain that an increase in the inclusions concentration $\phi$ tends to reduce the thickness of the convecting layer. Indeed, settling of the heavy inclusions occurs over a thicker layer (Fig. C1) while the melting front remains deeper (Fig. C2): both effects reduce the thickness of the convective layer. This in turn decreases the Rayleigh number of the unmolten region, reducing convection intensity and favouring the layering of light inclusions.

## 5.2.3 The crust's heterogeneity upon melting

Our Configurations 1 and 2 assume a relatively complex dependency of the rheology on melt fraction. However, such a switch in rheology from a temperature and strain-rate power-law relationship (eq. 6) to a melt fraction dependency (eq. 9) has become common practice in numerical geodynamics (e.g. Gerya et al. 2008; Perchuk et al. 2011; Schenker et al. 2012; Gerbault et al. 2018). Here we showed the key role of the crust’s heterogeneity in the formation of migmatite domes.

It would be interesting to carry a comparative study with other modelling approaches of partial melting continental crust. Such studies usually seek other purposes, such as mobile belt geodynamics in Archean and Proterozoic times (Sizova et al. 2015; Piccolo et al. 2019) or high pressure exhumation patterns in continental collision belts (Perchuk et al. 2011; Burov et al. 2014). Other theoretical models investigate multi-scale melt segregation processes (Räss et al. 2019; Schmeling et al. 2019; Culha et al. 2020; Schmeling et al. 2023) with advanced computational tools in order to deal with high resolution scales. It would be useful to investigate how these approaches also reproduce the segregation of partially molten felsic domes above a convective crust, in comparison to field estimates of dimensions and timing associated to migmatite terranes (e.g. Ledru et al. 2001; Riel et al. 2016; Vanderhaeghe et al. 2018). Such a comparison will help to confirm the requirements obtained in our study, regarding both the compositional heterogeneities and the boundary conditions.

## 5.2.4 The influence of the boundary conditions

The lateral boundaries were kept fixed in our models, consistent with the assumption that the orogenic system is in a transitional stage between continental convergence and post-orogenic collapse. Actually such a transitional stage may last several tens of millions
of years, as has been proposed in several ancient orogenic belts, based on geochronological data (Ganne et al. 2014; Laurent et al. 2018; Gerbault et al. 2018). This shows that the diapiric rise of partially molten rocks from the core of orogenic roots, for example driven by gravitational instabilities, could indeed be a standard process of the orogenic cycle. Nevertheless, lateral extension has also been shown to play a strong control on final exhumation through the more resistant upper crust (Whitney et al. 2004; Schenker et al. 2012). Several studies have tackled the competition between lateral 'tectonic' forces and buoyancy over the orogenic cycle, but this issue stands beyond the scope of this study. Top surface erosion and sedimentary processes were not accounted for here. However they obviously play a role in exhumation, by adding or removing a vertical force component that may not be uniform along the horizontal top surface, as has been shown by numerous studies before (e.g. Avouac & Burov 1996; Willett 1999; Burov et al. 2003; Gerbault et al. 2003). Hence these processes should be considered if a more precise model of the exhumation rates of specific domes is carried out.

In terms of numerical boundary conditions, periodic conditions (material and temperature) were imposed on the lateral boundaries; we showed their negligible impact on the development of convective cells (their characteristic width) in Louis-Napoléon et al. (2022).

Finally, commenting on our thermal boundary conditions at the base of the models; we made two strong assumptions: (1) a sudden heating from 600 to 1000 °C at the model initiation time and (2) a progressive cooling back to 600 °C from 5 Myr to 10 or 15 Myr on-wards. For approximation (1), we assumed that lithosphere slab-breakoff, delamination or flowing underlying magmas occurs below the crust within several million years, consistent with van Hunen & Allen (2011) and Duretz & Gerya (2013) (see also discussion in Garzanti et al. 2018). One must consider that although the initiation of slab break-off or delamination takes time, extra heat actually reaches the base of the orogenic crust only at its ultimate stage. Alternatively, basal heating may also occur via underplating of basaltic magma from the mantle, especially in subduction zone mountain belts. In any case this process can be considered to occur much faster than the subsequent cooling event (cf. approximation 2). Approximation (2) in turn, for example cooling assumed to occur over ca. 5 to 10 Myr, is consistent with Davies & von Blanckenburg (1995), but could actually occur over a much longer period (e.g. Ueda et al. 2012). Actually, application of cooling from below causes an inversion of the temperature profile such that heat tends to flow downwards from the crust into the mantle. This leads the base of the crust to become more viscous than the middle of the crust, which may potentially induce other gravity instabilities (not occurring here since the black inclusions remain denser than the BED). Alternative boundary conditions to simulate orogenic cooling would have been to reduce the thermal flow instead of temperature (e.g. Thompson & Connolly, 1995), or to stretch the bottom boundary, which requires to increase the model's size. Nonetheless, let us mention that a lithosphere's thermal state during a continental orogenic cycle is certainly not uniform neither in space nor time, and that lateral mass transfer within the crust (lateral flow) and below it (resuming slab advance as opposed to slab retreat) may lead to insulating layers at different levels of lithospheric columns.

Despite their imperfections, the thermal boundary conditions of our models are shown to strongly control the development and preservation of migmatite domes. Shorter heating time would tend to prevent convective motion, and longer heating time would inhibit diapirism. Hence they stand as critical conditions for dome

formation. We acknowledge that the global thermal evolution of an orogenic system should then better be tackled with larger scale models, that incorporate the entire lithosphere and asthenospheric domains and can then account for self-consistent lithospheric delamination and/or slab-break off processes (e.g. Schott et al. 2000; Ueda et al. 2012; Burov et al. 2014). However such large scale models also require to consider the orogens entire evolutionary history, requiring more precise thermal constraints from field data in order to validate specific geodynamic scenarios. They also remain submitted to limitations imposed by mesh resolution.

## 5.3 Insight on the conditions for the formation of migmatite domes worldwide

One may wonder how to transpose the controlling parameters displayed by our models to field observations. The answer is twofold, in terms of the observable spatial dimensions and time-scales first, and then in terms of the implications on the rheological properties and the boundary conditions of the system. First, domes of sizes 10–100 km, such as in the Archean Pilbara craton, the Palaeoproterozoic Eburnean and Svecofennian belts and even in the Palaeozoic Variscan belt (Velay area), likely formed in warmer crusts, consistent with a higher geothermal gradient (Brown & Johnson 2018; Vanderhaeghe et al. 2019) and/or the persistence of thick orogenic crust for several tens of Myr (Harley 2016; Laurent et al. 2018; Guergouz et al. 2018; Turlin et al. 2018). As of small scale domes, such as those of Thor-Odin in the Shuswap MCC or Naxos in the Aegean domain (Kruckenberg et al. 2011; Vanderhaeghe et al. 2018), the models consistently indicate that thermal convection and diapirism occurred over a relatively short time-window of 5–10 Myr. This timing indicates a relatively short thermal destabilization process, such as slab-break-off or lithospheric delamination followed by rapid thermal re-equilibration, prior to a drastic change in the lateral tectonic forces (which controls subsequent orogenic evolution, e.g. in the Aegean domain).

Second, our models confirm that an orogenic crust can convect and develop dome structures considering a reasonable range of density contrasts (100 kg m⁻³) and viscosities locally as low as 10¹⁶ Pa s. This mixing process (convection) followed by gravitationally driven segregation, occurs over several millions years, and may explain the magmas variabilities that are found in migmatite domes worldwide, with cyclic aggregation of molten material coming from below. A closer comparison between models and data requires further knowledge from the geological record: (i) measured field lineations can provide indications on the symmetry of the domes shapes, (ii) zircon dating and other geochronology data can provide a more precise timing of vertical (P,T) motions from the centre of the domes to their borders and (iii) mineralogical and fluid inclusions studies can help identify the (P, T, X) conditions associated with the ongoing metamorphic reactions, melting and crystallization, during the medium's deformation (from the domes core to their edges, and from deep to shallow depths). This in turn, can help to better identify the appropriate density and viscosity of elementary volumes representative of the physical scales at play (and link the grain scale and the crustal scale).

## 6 CONCLUSION

Our models show the general importance of considering the compositional heterogeneity of continental crust upon melting. Melting processes involve thermodynamic exchanges at multiple scales, and
of years, as has been proposed in several ancient orogenic belts, based on geochronological data (Ganne et al. 2014; Laurent et al. 2018; Gerbault et al. 2018). This shows that the diapiric rise of partially molten rocks from the core of orogenic roots, for example driven by gravitational instabilities, could indeed be a standard process of the orogenic cycle. Nevertheless, lateral extension has also been shown to play a strong control on final exhumation through the more resistant upper crust (Whitney et al. 2004; Schenker et al. 2012). Several studies have tackled the competition between lateral ‘tectonic’ forces and buoyancy over the orogenic cycle, but this issue stands beyond the scope of this study. Top surface erosion and sedimentary processes were not accounted for here. However they obviously play a role in exhumation, by adding or removing a vertical force component that may not be uniform along the horizontal top surface, as has been shown by numerous studies before (e.g. Avouac & Burov 1996; Willett 1999; Burov et al. 2003; Gerbault et al. 2003). Hence these processes should be considered if a more precise model of the exhumation rates of specific domes is carried out.

In terms of numerical boundary conditions, periodic conditions (material and temperature) were imposed on the lateral boundaries; we showed their negligible impact on the development of convective cells (their characteristic width) in Louis-Napoléon et al. (2022).

Finally, commenting on our thermal boundary conditions at the base of the models; we made two strong assumptions: (1) a sudden heating from 600 to 1000 °C at the model initiation time and (2) a progressive cooling back to 600 °C from 5 Myr to 10 or 15 Myr on-wards. For approximation (1), we assumed that lithosphere slab-breakoff, delamination or flowing underlying magmas occurs below the crust within several million years, consistent with van Hunen & Allen (2011) and Duretz & Gerya (2013) (see also discussion in Garzanti et al. 2018). One must consider that although the initiation of slab break-off or delamination takes time, extra heat actually reaches the base of the orogenic crust only at its ultimate stage. Alternatively, basal heating may also occur via underplating of basaltic magma from the mantle, especially in subduction zone mountain belts. In any case this process can be considered to occur much faster than the subsequent cooling event (cf. approximation 2). Approximation (2) in turn, for example cooling assumed to occur over ca. 5 to 10 Myr, is consistent with Davies & von Blanckenburg (1995), but could actually occur over a much longer period (e.g. Ueda et al. 2012). Actually, application of cooling from below causes an inversion of the temperature profile such that heat tends to flow downwards from the crust into the mantle. This leads the base of the crust to become more viscous than the middle of the crust, which may potentially induce other gravity instabilities (not occurring here since the black inclusions remain denser than the BED). Alternative boundary conditions to simulate orogenic cooling would have been to reduce the thermal flow instead of temperature (e.g. Thompson & Connolly, 1995), or to stretch the bottom boundary, which requires to increase the model's size. Nonetheless, let us mention that a lithosphere's thermal state during a continental orogenic cycle is certainly not uniform neither in space nor time, and that lateral mass transfer within the crust (lateral flow) and below it (resuming slab advance as opposed to slab retreat) may lead to insulating layers at different levels of lithospheric columns.

Despite their imperfections, the thermal boundary conditions of our models are shown to strongly control the development and preservation of migmatite domes. Shorter heating time would tend to prevent convective motion, and longer heating time would inhibit diapirism. Hence they stand as critical conditions for dome

formation. We acknowledge that the global thermal evolution of an orogenic system should then better be tackled with larger scale models, that incorporate the entire lithosphere and asthenospheric domains and can then account for self-consistent lithospheric delamination and/or slab-break off processes (e.g. Schott et al. 2000; Ueda et al. 2012; Burov et al. 2014). However such large scale models also require to consider the orogens entire evolutionary history, requiring more precise thermal constraints from field data in order to validate specific geodynamic scenarios. They also remain submitted to limitations imposed by mesh resolution.

## 5.3 Insight on the conditions for the formation of migmatite domes worldwide

One may wonder how to transpose the controlling parameters displayed by our models to field observations. The answer is twofold, in terms of the observable spatial dimensions and time-scales first, and then in terms of the implications on the rheological properties and the boundary conditions of the system. First, domes of sizes 10–100 km, such as in the Archean Pilbara craton, the Palaeoproterozoic Eburnean and Svecofennian belts and even in the Palaeozoic Variscan belt (Velay area), likely formed in warmer crusts, consistent with a higher geothermal gradient (Brown & Johnson 2018; Vanderhaeghe et al. 2019) and/or the persistence of thick orogenic crust for several tens of Myr (Harley 2016; Laurent et al. 2018; Guergouz et al. 2018; Turlin et al. 2018). As of small scale domes, such as those of Thor-Odin in the Shuswap MCC or Naxos in the Aegean domain (Kruckenberg et al. 2011; Vanderhaeghe et al. 2018), the models consistently indicate that thermal convection and diapirism occurred over a relatively short time-window of 5–10 Myr. This timing indicates a relatively short thermal destabilization process, such as slab-break-off or lithospheric delamination followed by rapid thermal re-equilibration, prior to a drastic change in the lateral tectonic forces (which controls subsequent orogenic evolution, e.g. in the Aegean domain).

Second, our models confirm that an orogenic crust can convect and develop dome structures considering a reasonable range of density contrasts (100 kg m⁻³) and viscosities locally as low as 10¹⁶ Pa s. This mixing process (convection) followed by gravitationally driven segregation, occurs over several millions years, and may explain the magmas variabilities that are found in migmatite domes worldwide, with cyclic aggregation of molten material coming from below. A closer comparison between models and data requires further knowledge from the geological record: (i) measured field lineations can provide indications on the symmetry of the domes shapes, (ii) zircon dating and other geochronology data can provide a more precise timing of vertical (P,T) motions from the centre of the domes to their borders and (iii) mineralogical and fluid inclusions studies can help identify the (P, T, X) conditions associated with the ongoing metamorphic reactions, melting and crystallization, during the medium's deformation (from the domes core to their edges, and from deep to shallow depths). This in turn, can help to better identify the appropriate density and viscosity of elementary volumes representative of the physical scales at play (and link the grain scale and the crustal scale).

## 6 CONCLUSION

Our models show the general importance of considering the compositional heterogeneity of continental crust upon melting. Melting processes involve thermodynamic exchanges at multiple scales, and
several geodynamic studies in the community now incorporate such multiscale mechanisms, via the implementation of a two-phase flow formalism. Here we approximated the subscale melt segregation processes by their first order properties, for example heterogeneity at the intermediate scale of several hundred metres. This allowed us to assess how gravitational dynamics together with the heating history can lead to the development and preservation of migmatite domes.

In our study part I (Louis-Napoléon et al. 2022), we showed that when convection occurs in partially molten crust, it tends to destroy any segregated dome made of the lighter products of partial melting. We had found an intermediate regime in which buoyant partially molten rocks can accumulate and form a horizontal ('stratified') layer at the base of the upper crust.

Here in this part II, we complemented this geodynamic context by exploring additional factors allowing for the formation and preservation of dome shaped structures in a convecting system. We achieved this by accounting for the following conditions:

(i) Basal heating for about 5–10 Myr prior to the return to standard basal temperature during another 5–10 Myr. A faster evolution of this thermal condition impedes the development of convection, while a slower evolution impedes the preservation of kilometric scale dome structures. Such timing is consistent with estimates of the typical timing of lithospheric delamination or break-off during orogenic evolution.

(ii) Non-reversible density of partially molten material when reaching its peak metamorphic grade facilitates the preservation of dome structures instead of layered structures (at the base of the upper crust).

(iii) In cases of a vigorous convective regime in the lower crust ('suspension' regime, Configuration 1), domes of sizes of about 10 km or more can form. This regime applies to a 'hot' felsic crust (Louis-Napoléon et al. 2022). Motion velocities reach several cm/yr.

(iv) In contrast when the system is set in the 'layering' regime (Louis-Napoléon et al. (2022), for example the bulk crust is more viscous (Configuration 2), smaller scale dome structures are preserved, of sizes rather close to 5 km. Moreover, these small domes have a more homogeneous composition made of buoyant material, in comparison to those of Configuration 1. In this regime, the system achieves moderate velocities (about 0.5 cm yr⁻¹ or less) that are consistent with present-day transitional tectonics velocities (switch from compressional to tensile tectonics).

Further work is still required to improve our understanding of the constraints on migmatite dome formation, and future models should account for (i) moving the lateral boundaries and considering erosion/sedimentation surface processes on more global time and space scales, (ii) small scale segregation processes that depend on compositional thermodynamics. Although our models reproduced both small (≤5) km and large (≥20 km) dome sizes, coherent with cooler orogenic settings in the Phanerozoic than in the Archean, an improved model fit to data requires further information on the detailed chronology and compositions associated with those terranes. Geological mapping and mineralogical compositions at different scales will help to better constrain and extrapolate properties known at the grain scale (transport of volatile species), to the metre (compaction and fracturing processes) and kilometric scales (mixed heterogeneous or homogeneous textures), and finally to the crustal scale. Precise geophysical imaging would also help visualize how these domes are actually anchored at depth, and where their heavier, deeper counter-parts have gone.

# SUPPORTING INFORMATION

Supplementary data are available at GJI online.

## suppl.data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

# ACKNOWLEDGMENTS

We thank the Institut National Polytechnique de Toulouse (INPT) which supported this work. Complementary financial support by INSU-Syster is acknowledged. This work was granted access to the HPC resources of CALMIP supercomputing centre under the allocation 2022-[1525]. Detailed and careful reviews by Harro Schmeling and Roberto Weinberg are warmly thanked, as they greatly helped clarifying and correcting many key points in this work.

# DATA AVAILABILITY

Our numerical models are performed with our own solver called *multiMeltInterFoamv2*, which is available here: https://gitlab.com/AurelieLN/MultiMeltInterFoam.git. Videos of the 3-D runs can be found here: https://www.youtube.com/ playlist?list=PLCdeVmWEr uk22N4GDS5mrEk0nNkJbTwFu

# REFERENCES

* Annen, C. & Sparks, R., 2002. Effects of repetitive emplacement of basaltic intrusions on thermal evolution and melt generation in the crust, *Earth planet. Sci. Lett.*, **203**(3–4), 937–955.

* Arzi, A.A., 1978. Critical phenomena in the rheology of partially melted rocks, *Tectonophysics*, **44**(1–4), 173–184.

* Avouac, J.-P. & Burov, E., 1996. Erosion as a driving mechanism of intracontinental mountain growth, *J. geophys. Res.*, **101**(B8), 17 747–17 769.

* Babeyko, A., Sobolev, S., Trumbull, R.B., O.O. & Lavier, L., 2002. Numerical models of crustal scale convection and partial melting beneath the Altiplano-Puna Plateau, *Earth planet. Sci. Lett.*, **199**(3), 373–388.

* Bachmann, O. & Huber, C., 2016. Silicic magma reservoirs in the Earth's crust, *Am. Mineral.*, **101**(11), 2377–2404.

* Balen, D. & Broska, I., 2011. Tourmaline nodules: products of devolatilization within the final evolutionary stage of granitic melt?, *Geol. Soc., Lond., Spec. Publ.*, **350**(1), 53–68.

* Berberović, E., van Hinsberg, N.P., Jakirlić, S., Roisman, I.V. & Tropea, C., 2009. Drop impact onto a liquid layer of finite thickness: dynamics of the cavity evolution, *Phys. Rev. E*, **79**, 036306, doi:10.1103/PhysRevE.79.036306.

* Bianco, M.J., Gerstoft, P., Olsen, K.B. & Lin, F.-C., 2019. High-resolution seismic tomography of Long Beach, CA using machine learning, *Sci. Rep.*, **9**(1), doi:10.1038/s41598-019-50381-z.

* Billerot, A., Duchene, S., Vanderhaeghe, O. & de Sigoyer, J., 2017. Gneiss domes of the Danba metamorphic complex, Songpan Ganze, Eastern Tibet, *J. Asian Earth Sci.*, **140**, 48–74.

* Brown, M., 1995. Mechanisms and consequence of melt segregation from crustal protolith, *J. geophys. Res.*, **100**, 15 655–15 679.

* Brown, M., 2010. Melting of the continental crust during orogenesis: the thermal, rheological, and compositional consequences of melt transport from lower to upper continental crust, *Can. J. Earth Sci.*, **47**(5), 655–694.

* Brown, M. & Johnson, T., 2018. Secular change in metamorphism and the onset of global plate tectonics, *Am. Mineral.: J. Earth planet. Mater.*, **103**(2), 181–196.
several geodynamic studies in the community now incorporate such multiscale mechanisms, via the implementation of a two-phase flow formalism. Here we approximated the subscale melt segregation processes by their first order properties, for example heterogeneity at the intermediate scale of several hundred metres. This allowed us to assess how gravitational dynamics together with the heating history can lead to the development and preservation of migmatite domes.

In our study part I (Louis-Napoléon et al. 2022), we showed that when convection occurs in partially molten crust, it tends to destroy any segregated dome made of the lighter products of partial melting. We had found an intermediate regime in which buoyant partially molten rocks can accumulate and form a horizontal ('stratified') layer at the base of the upper crust.

Here in this part II, we complemented this geodynamic context by exploring additional factors allowing for the formation and preservation of dome shaped structures in a convecting system. We achieved this by accounting for the following conditions:

(i) Basal heating for about 5–10 Myr prior to the return to standard basal temperature during another 5–10 Myr. A faster evolution of this thermal condition impedes the development of convection, while a slower evolution impedes the preservation of kilometric scale dome structures. Such timing is consistent with estimates of the typical timing of lithospheric delamination or break-off during orogenic evolution.

(ii) Non-reversible density of partially molten material when reaching its peak metamorphic grade facilitates the preservation of dome structures instead of layered structures (at the base of the upper crust).

(iii) In cases of a vigorous convective regime in the lower crust ('suspension' regime, Configuration 1), domes of sizes of about 10 km or more can form. This regime applies to a 'hot' felsic crust (Louis-Napoléon et al. 2022). Motion velocities reach several cm/yr.

(iv) In contrast when the system is set in the 'layering' regime (Louis-Napoléon et al. (2022), for example the bulk crust is more viscous (Configuration 2), smaller scale dome structures are preserved, of sizes rather close to 5 km. Moreover, these small domes have a more homogeneous composition made of buoyant material, in comparison to those of Configuration 1. In this regime, the system achieves moderate velocities (about 0.5 cm yr⁻¹ or less) that are consistent with present-day transitional tectonics velocities (switch from compressional to tensile tectonics).

Further work is still required to improve our understanding of the constraints on migmatite dome formation, and future models should account for (i) moving the lateral boundaries and considering erosion/sedimentation surface processes on more global time and space scales, (ii) small scale segregation processes that depend on compositional thermodynamics. Although our models reproduced both small (≤5) km and large (≥20 km) dome sizes, coherent with cooler orogenic settings in the Phanerozoic than in the Archean, an improved model fit to data requires further information on the detailed chronology and compositions associated with those terranes. Geological mapping and mineralogical compositions at different scales will help to better constrain and extrapolate properties known at the grain scale (transport of volatile species), to the metre (compaction and fracturing processes) and kilometric scales (mixed heterogeneous or homogeneous textures), and finally to the crustal scale. Precise geophysical imaging would also help visualize how these domes are actually anchored at depth, and where their heavier, deeper counter-parts have gone.

# SUPPORTING INFORMATION

Supplementary data are available at GJI online.

## suppl.data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

# ACKNOWLEDGMENTS

We thank the Institut National Polytechnique de Toulouse (INPT) which supported this work. Complementary financial support by INSU-Syster is acknowledged. This work was granted access to the HPC resources of CALMIP supercomputing centre under the allocation 2022-[1525]. Detailed and careful reviews by Harro Schmeling and Roberto Weinberg are warmly thanked, as they greatly helped clarifying and correcting many key points in this work.

# DATA AVAILABILITY

Our numerical models are performed with our own solver called *multiMeltInterFoamv2*, which is available here: https://gitlab.com/AurelieLN/MultiMeltInterFoam.git. Videos of the 3-D runs can be found here: https://www.youtube.com/ playlist?list=PLCdeVmWEruk22N4GDS5mrEk0nNkJbTwFu

# REFERENCES

* Annen, C. & Sparks, R., 2002. Effects of repetitive emplacement of basaltic intrusions on thermal evolution and melt generation in the crust, *Earth planet. Sci. Lett.*, **203**(3–4), 937–955.

* Arzi, A.A., 1978. Critical phenomena in the rheology of partially melted rocks, *Tectonophysics*, **44**(1–4), 173–184.

* Avouac, J.-P. & Burov, E., 1996. Erosion as a driving mechanism of intracontinental mountain growth, *J. geophys. Res.*, **101**(B8), 17 747–17 769.

* Babeyko, A., Sobolev, S., Trumbull, R.B., O.O. & Lavier, L., 2002. Numerical models of crustal scale convection and partial melting beneath the Altiplano-Puna Plateau, *Earth planet. Sci. Lett.*, **199**(3), 373–388.

* Bachmann, O. & Huber, C., 2016. Silicic magma reservoirs in the Earth's crust, *Am. Mineral.*, **101**(11), 2377–2404.

* Balen, D. & Broska, I., 2011. Tourmaline nodules: products of devolatilization within the final evolutionary stage of granitic melt?, *Geol. Soc., Lond., Spec. Publ.*, **350**(1), 53–68.

* Berberović, E., van Hinsberg, N.P., Jakirlić, S., Roisman, I.V. & Tropea, C., 2009. Drop impact onto a liquid layer of finite thickness: dynamics of the cavity evolution, *Phys. Rev. E*, **79**, 036306, doi:10.1103/PhysRevE.79.036306.

* Bianco, M.J., Gerstoft, P., Olsen, K.B. & Lin, F.-C., 2019. High-resolution seismic tomography of Long Beach, CA using machine learning, *Sci. Rep.*, **9**(1), doi:10.1038/s41598-019-50381-z.

* Billerot, A., Duchene, S., Vanderhaeghe, O. & de Sigoyer, J., 2017. Gneiss domes of the Danba metamorphic complex, Songpan Ganze, Eastern Tibet, *J. Asian Earth Sci.*, **140**, 48–74.

* Brown, M., 1995. Mechanisms and consequence of melt segregation from crustal protolith, *J. geophys. Res.*, **100**, 15 655–15 679.

* Brown, M., 2010. Melting of the continental crust during orogenesis: the thermal, rheological, and compositional consequences of melt transport from lower to upper continental crust, *Can. J. Earth Sci.*, **47**(5), 655–694.

* Brown, M. & Johnson, T., 2018. Secular change in metamorphism and the onset of global plate tectonics, *Am. Mineral.: J. Earth planet. Mater.*, **103**(2), 181–196.
Brown, M., Averkin, Y.A., McLellan, E.L. & Sawyer, E.W., 1995. Melt segregation in migmatites, *J. geophys. Res.*, **100**(B8), 15 655–15 679.

Brun, J.-P., 1983. L'origine des domes gneissiques; modeles et tests, *Bulletin de la Société géologique de France*, **7**(2), 219–228.

Buck, W.R., 1991. Modes of continental lithospheric extension, *J. geophys. Res.*, **96**(B12), 20 161–20 178.

Burg, J.-P. & Vanderhaeghe, O., 1993. Structures and way-up criteria in migmatites, with application to the Velay dome (French Massif Central), *J. Struct. Geol.*, **15**(11), 1293–1301.

Burg, J.-P., Kaus, B.J. & Podladchikov, Y.Y., 2004. Dome structures in collision orogens: mechanical investigation of the gravity/compression interplay, in *Gneiss Domes in Orogeny*, Vol. **380**, GSA Special Papers, pp. 47–66, eds Whitney, D.L., Teyssier, C. & Siddoway, C.S., Geological Society of America.

Burnham, C.W., 1979. The importance of volatile constituents, in *The Evolution of the Igneous Rocks*, pp. 439–482, Yoder, H. S., Princeton Univ. Press.

Burov, E., Jaupart, C. & Guillou-Frottier, L., 2003. Ascent and emplacement of buoyant magma bodies in brittle-ductile upper crust, *J. geophys. Res.*, **108**(B4), 2177, doi:10.1029/2002JB001904.

Burov, E. et al., 2014. Rheological and geodynamic controls on the mechanisms of subduction and HP/UHP exhumation of crustal rocks during continental collision: insights from numerical models, *Tectonophysics*, **631**, 212–250.

Caldwell, W.B., Klemperer, S.L., Rai, S.S. & Lawrence, J.F., 2009. Partial melt in the upper-middle crust of the Northwest Himalaya revealed by Rayleigh wave dispersion, *Tectonophysics*, **477**(1–2), 58–65.

Cao, W., Kaus, B. & Paterson, S., 2016. Intrusion of granitic magma into the continental crust facilitated by magma pulsing and dike-diapir interactions: numerical simulations, *Tectonics*, **35**(6), 1575–1594.

Cashman, K.V., Sparks, R. S.J. & Blundy, J.D., 2017. Vertically extensive and unstable magmatic systems: a unified view of igneous processes, *Science*, **355**(6331), doi:10.1126/science.aag3055.

Chen, J. et al., 2018. Melting conditions in the modern tibetan crust since the miocene, *Nat. Commun.*, **9**(1), 1–13.

Chen, Y. & Morgan, W.J., 1990. A nonlinear rheology model for mid-ocean ridge axis topography, *J. geophys. Res.*, **95**(B11), 17 583–17 604.

Chopin, F. et al., 2020. The Vaasa Migmatitic Complex (Svecofennian Orogen, Finland): Buildup of a LP-HT Dome During Nuna Assembly, *Tectonics*, **39**(3), doi:10.1029/2019TC005583.

Clos, F., Weinberg, R.F., Zibra, I. & Fenwick, M.J., 2019. Archean diapirism recorded by vertical sheath folds in the core of the yalgoo dome, yilgarn craton, *Precambrian Res.*, **320**, 391–402.

Collins, W.J., Van Kranendonk, M.J. & Teyssier, C., 1998. Partial convective overturn of Archaean crust in the east Pilbara Craton, western Australia: driving mechanisms and tectonic implications, *J. Struct. Geol.*, **20**(9–10), 1405–1424.

Collins, W., 1989. Polydiapirism of the Archean Mount Edgar Batholith, Pilbara Block, Western Australia, *Precambrian Res.*, **43**(1–2), 41–62.

Connolly, J. & Podladchikov, Y.Y., 1998. Compaction-driven fluid flow in viscoelastic rock, *Geodin. Acta*, **11**(2–3), 55–84.

Cruden, A., Koyi, H. & Schmeling, H., 1995. Diapiric basal entrainment of mafic into felsic magma, *Earth planet. Sci. Lett.*, **131**(3–4), 321–340.

Culha, C., Suckale, J., Keller, T. & Qin, Z., 2020. Crystal fractionation by crystal-driven convection, *Geophys. Res. Lett.*, **47**(4), e2019GL086784, doi:10.1029/2019GL086784.

Davies, J.H. & von Blanckenburg, F., 1995. Slab breakoff: a model of lithosphere detachment and its test in the magmatism and deformation of collisional orogens, *Earth planet. Sci. Lett.*, **129**(1–4), 85–102.

de Saint Blanquat, M., Horsman, E., Habert, G., Morgan, S., Vanderhaeghe, O., Law, R. & Tikoff, B., 2011. Multiscale magmatic cyclicity, duration of pluton construction, and the paradoxical relationship between tectonism and plutonism in continental arcs, *Tectonophysics*, **500**(1–4), 20–33.

Deubelbeiss, Y. & Kaus, B., 2008. Comparison of Eulerian and Lagrangian numerical techniques for the stokes equations in the presence of strongly varying viscosity, *Phys. Earth planet. Inter.*, **171**(1), 92–111.

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**(5), 301–321.

Duretz, T. & Gerya, T., 2013. Slab detachment during continental collision: influence of crustal rheology and interaction with lithospheric delamination, *Tectonophysics*, **602**, 124–140.

Edmonds, M., Cashman, K.V., Holness, M. & Jackson, M., 2019. Architecture and dynamics of magma reservoirs, *Phil. Trans. R. Soc.*, **A**, **377**, doi:10.1098/rsta.2018.0298.

Eskola, P.E., 1948. The problem of mantled gneiss domes, *Quart. J. Geol. Soc.*, **104**(1–4), 461–476.

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.

Garzanti, E., Radeff, G. & Malusà, M.G., 2018. Slab breakoff: a critical appraisal of a geological theory as applied in space and time, *Earth-Sci. Rev.*, **177**, 303–319.

Gerbault, M., Henrys, S. & Davey, F., 2003. Numerical models of lithospheric deformation forming the Southern Alps of New Zealand, *J. geophys. Res.*, **108**(B7), doi:10.1029/2001JB001716.

Gerbault, M., Schneider, J., Reverso-Peila, A. & Corsini, M., 2018. Crustal exhumation during ongoing compression in the Variscan Maures-Tanneron Massif, France-Geological and thermo-mechanical aspects, *Tectonophysics*, **746**, 439–458.

Gerya, T., Perchuk, L. & Burg, J.-P., 2008. Transient hot channels: perpetrating and regurgitating ultrahigh-pressure, high-temperature crust-mantle associations in collision belts, *Lithos*, **103**(1–2), 236–256.

Glazner, A.F., 2019. The ascent of water-rich magma and decompression heating: a thermodynamic analysis, *Am. Mineral.: J. Earth planet. Mater.*, **104**(6), 890–896.

Guergouz, C., Martin, L., Vanderhaeghe, O., Thébaud, N. & Fiorentini, M., 2018. Zircon and monazite petrochronologic record of prolonged amphibolite to granulite facies metamorphism in the Ivrea-Verbano and Strona-Ceneri Zones, NW Italy, *Lithos*, **308**, 1–18.

Harada, S., Mitsui, T. & Sato, K., 2012. Particle-like and fluid-like settling of a stratified suspension, *Eur. Phys. J. E*, **35**(1), 1–6.

Harley, S.L., 2016. A matter of time: the importance of the duration of UHT metamorphism, *J. Mineral. Petrol. Sci.*, **111**(2), 50–72.

Holloway, J.R., 1976. Fluids in the evolution of granitic magmas: consequences of finite CO₂ solubility, *Bull. geol. Soc. Am.*, **87**(10), 1513–1518.

Jasak, H., Weller, H. & Gosman, A., 1999. High resolution nvd differencing scheme for arbitrarily unstructured meshes, *Int. J. Numer. Methods Fluids*, **31**(2), 431–449.

Kisters, A., Ward, R., Anthonissen, C. & Vietze, M., 2009. Melt segregation and far-field melt transfer in the mid-crust, *J. Geol. Soc.*, **166**(5), 905–918.

Krieger, I.M. & Dougherty, T.J., 1959. A mechanism for non-Newtonian flow in suspensions of rigid spheres, *Trans. Soc. Rheol.*, **3**(1), 137–152.

Kruckenberg, S., Vanderhaeghe, O., Ferré, E., Teyssier, C. & Whitney, D., 2011. Flow of partially molten crust and the internal dynamics of a migmatite dome, Naxos, Greece, *Tectonics*, **30**(3), TC3001, doi:10.1029/2010TC002751.

Kruckenberg, S.C., Whitney, D.L., Teyssier, C., Fanning, C.M. & Dunlap, W.J., 2008. Paleocene-eocene migmatite crystallization, extension, and exhumation in the hinterland of the northern cordillera: Okanogan dome, washington, usa, *Bull. geol. Soc. Am.*, **120**(7–8), 912–929.

Laurent, A.T., Bingen, B., Duchene, S., Whitehouse, M.J., Seydoux-Guillaume, A.-m. & Bosse, V., 2018. Decoding a protracted zircon geochronological record in ultrahigh temperature granulite, and persistence of partial melting in the crust, Rogaland, Norway, *Contrib. Mineral. Petrol.*, **173**(4), 1–25.

Lavorel, G. & Le Bars, M., 2009. Sedimentation of particles in a vigorously convecting fluid, *Phys. Rev. E*, **80**(4), doi:10.1103/PhysRevE.80.046324.

Ledru, P., Courrioux, G., Dallain, C., Lardeaux, J., Montel, J., Vanderhaeghe, O. & Vitel, G., 2001. The Velay dome (French Massif Central): melt generation and granite emplacement during orogenic evolution, *Tectonophysics*, **342**(3–4), 207–237.
Brown, M., Averkin, Y.A., McLellan, E.L. & Sawyer, E.W., 1995. Melt segregation in migmatites, *J. geophys. Res.*, **100**(B8), 15 655–15 679.

Brun, J.-P., 1983. L'origine des domes gneissiques; modeles et tests, *Bulletin de la Société géologique de France*, **7**(2), 219–228.

Buck, W.R., 1991. Modes of continental lithospheric extension, *J. geophys. Res.*, **96**(B12), 20 161–20 178.

Burg, J.-P. & Vanderhaeghe, O., 1993. Structures and way-up criteria in migmatites, with application to the Velay dome (French Massif Central), *J. Struct. Geol.*, **15**(11), 1293–1301.

Burg, J.-P., Kaus, B.J. & Podladchikov, Y.Y., 2004. Dome structures in collision orogens: mechanical investigation of the gravity/compression interplay, in *Gneiss Domes in Orogeny*, Vol. **380**, GSA Special Papers, pp. 47–66, eds Whitney, D.L., Teyssier, C. & Siddoway, C.S., Geological Society of America.

Burnham, C.W., 1979. The importance of volatile constituents, in *The Evolution of the Igneous Rocks*, pp. 439–482, Yoder, H. S., Princeton Univ. Press.

Burov, E., Jaupart, C. & Guillou-Frottier, L., 2003. Ascent and emplacement of buoyant magma bodies in brittle-ductile upper crust, *J. geophys. Res.*, **108**(B4), 2177, doi:10.1029/2002JB001904.

Burov, E. et al., 2014. Rheological and geodynamic controls on the mechanisms of subduction and HP/UHP exhumation of crustal rocks during continental collision: insights from numerical models, *Tectonophysics*, **631**, 212–250.

Caldwell, W.B., Klemperer, S.L., Rai, S.S. & Lawrence, J.F., 2009. Partial melt in the upper-middle crust of the Northwest Himalaya revealed by Rayleigh wave dispersion, *Tectonophysics*, **477**(1–2), 58–65.

Cao, W., Kaus, B. & Paterson, S., 2016. Intrusion of granitic magma into the continental crust facilitated by magma pulsing and dike-diapir interactions: numerical simulations, *Tectonics*, **35**(6), 1575–1594.

Cashman, K.V., Sparks, R. S.J. & Blundy, J.D., 2017. Vertically extensive and unstable magmatic systems: a unified view of igneous processes, *Science*, **355**(6331), doi:10.1126/science.aag3055.

Chen, J. et al., 2018. Melting conditions in the modern tibetan crust since the miocene, *Nat. Commun.*, **9**(1), 1–13.

Chen, Y. & Morgan, W.J., 1990. A nonlinear rheology model for mid-ocean ridge axis topography, *J. geophys. Res.*, **95**(B11), 17 583–17 604.

Chopin, F. et al., 2020. The Vaasa Migmatitic Complex (Svecofennian Orogen, Finland): Buildup of a LP-HT Dome During Nuna Assembly, *Tectonics*, **39**(3), doi:10.1029/2019TC005583.

Clos, F., Weinberg, R.F., Zibra, I. & Fenwick, M.J., 2019. Archean diapirism recorded by vertical sheath folds in the core of the yalgoo dome, yilgarn craton, *Precambrian Res.*, **320**, 391–402.

Collins, W.J., Van Kranendonk, M.J. & Teyssier, C., 1998. Partial convective overturn of Archaean crust in the east Pilbara Craton, western Australia: driving mechanisms and tectonic implications, *J. Struct. Geol.*, **20**(9–10), 1405–1424.

Collins, W., 1989. Polydiapirism of the Archean Mount Edgar Batholith, Pilbara Block, Western Australia, *Precambrian Res.*, **43**(1–2), 41–62.

Connolly, J. & Podladchikov, Y.Y., 1998. Compaction-driven fluid flow in viscoelastic rock, *Geodin. Acta*, **11**(2–3), 55–84.

Cruden, A., Koyi, H. & Schmeling, H., 1995. Diapiric basal entrainment of mafic into felsic magma, *Earth planet. Sci. Lett.*, **131**(3–4), 321–340.

Culha, C., Suckale, J., Keller, T. & Qin, Z., 2020. Crystal fractionation by crystal-driven convection, *Geophys. Res. Lett.*, **47**(4), e2019GL086784, doi:10.1029/2019GL086784.

Davies, J.H. & von Blanckenburg, F., 1995. Slab breakoff: a model of lithosphere detachment and its test in the magmatism and deformation of collisional orogens, *Earth planet. Sci. Lett.*, **129**(1–4), 85–102.

de Saint Blanquat, M., Horsman, E., Habert, G., Morgan, S., Vanderhaeghe, O., Law, R. & Tikoff, B., 2011. Multiscale magmatic cyclicity, duration of pluton construction, and the paradoxical relationship between tectonism and plutonism in continental arcs, *Tectonophysics*, **500**(1–4), 20–33.

Deubelbeiss, Y. & Kaus, B., 2008. Comparison of Eulerian and Lagrangian numerical techniques for the stokes equations in the presence of strongly varying viscosity, *Phys. Earth planet. Inter.*, **171**(1), 92–111.

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**(5), 301–321.

Duretz, T. & Gerya, T., 2013. Slab detachment during continental collision: influence of crustal rheology and interaction with lithospheric delamination, *Tectonophysics*, **602**, 124–140.

Edmonds, M., Cashman, K.V., Holness, M. & Jackson, M., 2019. Architecture and dynamics of magma reservoirs, *Phil. Trans. R. Soc.*, **A**, **377**, doi:10.1098/rsta.2018.0298.

Eskola, P.E., 1948. The problem of mantled gneiss domes, *Quart. J. Geol. Soc.*, **104**(1–4), 461–476.

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.

Garzanti, E., Radeff, G. & Malusà, M.G., 2018. Slab breakoff: a critical appraisal of a geological theory as applied in space and time, *Earth-Sci. Rev.*, **177**, 303–319.

Gerbault, M., Henrys, S. & Davey, F., 2003. Numerical models of lithospheric deformation forming the Southern Alps of New Zealand, *J. geophys. Res.*, **108**(B7), doi:10.1029/2001JB001716.

Gerbault, M., Schneider, J., Reverso-Peila, A. & Corsini, M., 2018. Crustal exhumation during ongoing compression in the Variscan Maures-Tanneron Massif, France-Geological and thermo-mechanical aspects, *Tectonophysics*, **746**, 439–458.

Gerya, T., Perchuk, L. & Burg, J.-P., 2008. Transient hot channels: perpetrating and regurgitating ultrahigh-pressure, high-temperature crust-mantle associations in collision belts, *Lithos*, **103**(1–2), 236–256.

Glazner, A.F., 2019. The ascent of water-rich magma and decompression heating: a thermodynamic analysis, *Am. Mineral.: J. Earth planet. Mater.*, **104**(6), 890–896.

Guergouz, C., Martin, L., Vanderhaeghe, O., Thébaud, N. & Fiorentini, M., 2018. Zircon and monazite petrochronologic record of prolonged amphibolite to granulite facies metamorphism in the Ivrea-Verbano and Strona-Ceneri Zones, NW Italy, *Lithos*, **308**, 1–18.

Harada, S., Mitsui, T. & Sato, K., 2012. Particle-like and fluid-like settling of a stratified suspension, *Eur. Phys. J. E*, **35**(1), 1–6.

Harley, S.L., 2016. A matter of time: the importance of the duration of UHT metamorphism, *J. Mineral. Petrol. Sci.*, **111**(2), 50–72.

Holloway, J.R., 1976. Fluids in the evolution of granitic magmas: consequences of finite CO₂ solubility, *Bull. geol. Soc. Am.*, **87**(10), 1513–1518.

Jasak, H., Weller, H. & Gosman, A., 1999. High resolution nvd differencing scheme for arbitrarily unstructured meshes, *Int. J. Numer. Methods Fluids*, **31**(2), 431–449.

Kisters, A., Ward, R., Anthonissen, C. & Vietze, M., 2009. Melt segregation and far-field melt transfer in the mid-crust, *J. Geol. Soc.*, **166**(5), 905–918.

Krieger, I.M. & Dougherty, T.J., 1959. A mechanism for non-Newtonian flow in suspensions of rigid spheres, *Trans. Soc. Rheol.*, **3**(1), 137–152.

Kruckenberg, S., Vanderhaeghe, O., Ferré, E., Teyssier, C. & Whitney, D., 2011. Flow of partially molten crust and the internal dynamics of a migmatite dome, Naxos, Greece, *Tectonics*, **30**(3), TC3001, doi:10.1029/2010TC002751.

Kruckenberg, S.C., Whitney, D.L., Teyssier, C., Fanning, C.M. & Dunlap, W.J., 2008. Paleocene-eocene migmatite crystallization, extension, and exhumation in the hinterland of the northern cordillera: Okanogan dome, washington, usa, *Bull. geol. Soc. Am.*, **120**(7–8), 912–929.

Laurent, A.T., Bingen, B., Duchene, S., Whitehouse, M.J., Seydoux-Guillaume, A.-m. & Bosse, V., 2018. Decoding a protracted zircon geochronological record in ultrahigh temperature granulite, and persistence of partial melting in the crust, Rogaland, Norway, *Contrib. Mineral. Petrol.*, **173**(4), 1–25.

Lavorel, G. & Le Bars, M., 2009. Sedimentation of particles in a vigorously convecting fluid, *Phys. Rev. E*, **80**(4), doi:10.1103/PhysRevE.80.046324.

Ledru, P., Courrioux, G., Dallain, C., Lardeaux, J., Montel, J., Vanderhaeghe, O. & Vitel, G., 2001. The Velay dome (French Massif Central): melt generation and granite emplacement during orogenic evolution, *Tectonophysics*, **342**(3–4), 207–237.
Lehujeur, M. & Chevrot, S., 2020. Eikonal tomography using coherent surface waves extracted from ambient noise by iterative matched filtering application to the large-n Maupasacq Array, *J. geophys. Res.*, **125**(6), e2020JB019363, doi:10.1029/2020JB019363.

Leitch, A. & Weinberg, R.F., 2002. Modelling granite migration by mesoscale pervasive flow, *Earth planet. Sci. Lett.*, **200**(1–2), 131–146.

Louis-Napoléon, A., 2020. Modélisation d'instabilités gravitaires au sein de la croûte partielle fondue par une méthode de Volume-Of-Fluid, *PhD thesis*, Toulouse, INPT, France.

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**(1), 474–506.

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

Parmigiani, A., Huber, C. & Bachmann, O., 2014. Mush microphysics and the reactivation of crystal-rich magma reservoirs, *J. geophys. Res.*, **119**(8), 6308–6322.

Patočka, V., Calzavarini, E. & Tosi, N., 2020. Settling of inertial particles in turbulent Rayleigh-Bénard convection, *Phys. Rev. Fluids*, **5**(11), doi:10.1103/PhysRevFluids.5.114304.

Perchuk, L.L., Gerya, T.V., van Reenen, D., Kramers, J. & McCourt, S., 2011. Formation and evolution of precambrian granulite terranes: a gravitational redistribution model, *Geol. Soc. Am. Memoirs*, **207**, 289–310.

Petford, N., 1996. Dykes or diapirs?, *Earth Environ. Sci. Trans. R. Soc. Edinburgh*. **87**(1–2), 105–114.

Piccolo, A., Palin, R.M., Kaus, B.J. & White, R.W., 2019. Generation of Earth's early continents from a relatively cool Archean mantle, *Geochem. Geophys. Geosyst.*, **20**(4), 1679–1697.

Pinkerton, H. & Stevenson, R., 1992. Methods of determining the rheological properties of magmas at sub-liquidus temperatures, *J. Volc. Geotherm. Res.*, **53**(1), 47–66.

Platt, J., 1980. Archaean greenstone belts: a structural test of tectonic hypothesis, *Tectonophysics*, **65**(1–2), 127–150.

Ramberg, H., 1981a. The role of gravity in orogenic belts, *Geol. Soc., Lond., Spec. Publ.*, **9**(1), 125–140.

Ramberg, H., 1981b. *Gravity, Deformation, and the Earth's Crust: In Theory, Experiments, and Geological Application*, 2nd edn, Academic Press, 452pp.

Ramsay, J.G., 1967. *Folding and Fracturing of Rocks*, Mc Graw Hill Book Company, 568pp.

Ranalli, G., 1995. *Rheology of the Earth*, Springer Science & Business Media.

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

Rey, P.F., Teyssier, C., Kruckenberg, S.C. & Whitney, D.L., 2011. Viscous collision in channel explains double domes in metamorphic core complexes, *Geol.*, **39**(4), 387–390.

Riel, N., Mercier, J. & Weinberg, R., 2016. Convection in a partially molten metasedimentary crust? Insights from the El Oro complex (Ecuador), *Geology*, **44**(1), 31–34.

Roberts, N.M. & Tikoff, B., 2021. Internal structure of the Paleoarchean Mt Edgar Dome, Pilbara Craton, Western Australia, *Precambrian Res.*, **358**, doi:10.1016/j.precamres.2021.106163.

Roscoe, R., 1952. The viscosity of suspensions of rigid spheres, *Br. J. Appl. Phys.*, **3**(8), doi:10.1088/0508-3443/3/8/306.

Rosenberg, C. & Handy, M., 2005. Experimental deformation of partially melted granite revisited: implications for the continental crust, *J. Metamor. Geol.*, **23**(1), 19–28.

Rutter, E., Brodie, K. & Irving, D., 2006. Flow of synthetic, wet, partially molten ‘granite’ under undrained conditions: an experimental study, *J. geophys. Res.*, **111**(B6), doi:10.1029/2005JB004257.

Sawyer, E., 1994. Melt segregation in the continental crust, *Geology*, **22**(11), 1019–1022.

Sawyer, E., Dombrowski, C. & Collins, W., 1999. Movement of melt during synchronous regional deformation and granulite-facies anatexis, an example from the Wuluma Hills, Central Australia, *Geol. Soc. Lond., Spec. Pub.*, **168**(1), 221–237.

Schenker, Filippo, L., 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.

Schilling, F.R. et al., 2006. Partial melting in the Central Andean Crust: a review of geophysical, petrophysical, and petrologic evidence, in *The Andes*, pp. 459–474, Oncken, O., et al., Springer.

Schmeling, H., Kruse, J.P. & Richard, G., 2012. Effective shear and bulk viscosity of partially molten rock based on elastic moduli theory of a fluid filled poroelastic medium, *Geophys. J. Int.*, **190**(3), 1571–1578.

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**(1), 422–450.

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**(5), e2023GC010860, doi:10.1029/2023GC010860.

Schott, B., Yuen, D.A. & Schmeling, H., 2000. The significance of shear heating in continental delamination, *Phys. Earth planet. Inter.*, **118**(3–4), 273–290.

Sizova, E., Gerya, T., Stüwe, K. & Brown, M., 2015. Generation of felsic crust in the Archean: a geodynamic modeling perspective, *Precambrian Res.*, **271**, 198–224.

Solar, G.S. & Brown, M., 2001. Deformation partitioning during transpression in response to Early Devonian oblique convergence, Northern Appalachian Orogen, USA, *J. Struct. Geol.*, **23**(6–7), 1043–1065.

Spang, A., Baumann, T.S. & Kaus, B.J., 2021. A multiphysics approach to constrain the dynamics of the Altiplano-Puna magmatic system, *J. geophys. Res.*, **126**(7), e2021JB021725, doi:10.1029/2021JB021725.

Thompson, A.B. & Connolly, J.A., 1995. Melting of the continental crust: some thermal and petrological constraints on anatexis in continental collision zones and other tectonic settings, *J. geophys. Res.*, **100**(B8), 1565–1579.

Toé, W., Vanderhaeghe, O., André-Mayer, A.-S., Feybesse, J.-L. & Milési, J.-P., 2013. From migmatites to granites in the Pan-African Damara orogenic belt, *J. Afr. Earth Sci.*, **85**, 62–74.

Turlin, F., Deruy, C., Eglinger, A., Vanderhaeghe, O., André-Mayer, A.-S., Poujol, M., Moukhsil, A. & Solgadi, F., 2018. A 70 ma record of suprasolidus conditions in the large, hot, long-duration Grenville Orogen, *Terra Nova*, **30**, 233–243.

Ueda, K., Gerya, T.V. & Burg, J., 2012. Delamination in collisional orogens: thermomechanical modeling, *J. geophys. Res.*, **117**(B8), doi:10.1029/2012JB009144.

Ulvrová, M., Labrosse, S., Cottice, N., Ráback, P. & Tackley, P., 2012. Numerical modelling of convection interacting with a melting and solidification front: application to the thermal evolution of the basal magma ocean, *Phys. Earth planet. Inter.*, **206**, 51–66.

Van der Molen, I. & Paterson, M., 1979. Experimental deformation of partially-melted granite, *Contrib. Mineral. Petrol.*, **70**(3), 299–318.

van Hunen, J. & Allen, M.B., 2011. Continental collision and slab break-off: a comparison of 3-D numerical models with observations, *Earth planet. Sci. Lett.*, **302**(1–2), 27–37.

Van Kranendonk, M.J., 2011. Cool greenstone drips and the role of partial convective overturn in Barberton Greenstone belt evolution, *J. Afr. Earth Sci.*, **60**(5), 346–352.

Van Kranendonk, M.J., Collins, W., Hickman, A. & Pawley, M.J., 2004. Critical tests of vertical vs. horizontal tectonic models for the Archaean east Pilbara Granite-Greenstone Terrane, Pilbara Craton, Western Australia, *Precambrian Res.*, **131**(3–4), 173–211.

Van Kranendonk, M.J., Hugh Smithies, R., Hickman, A.H. & Champion, D., 2007. secular tectonic evolution of archean continental crust: interplay