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
between horizontal and vertical processes in the formation of the Pilbara Craton, Australia, *Terra Nova*, **19**(1), 1–38.

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

Vanderhaeghe, O. & Teyssier, C., 2001. Partial melting and flow of orogens, *Tectonophysics*, **342**(3–4), 451–472.

Vanderhaeghe, O., Teyssier, C. & Wysoczanski, R., 1999. Structural and geochronological constraints on the role of partial melting during the formation of the Shuswap metamorphic core complex at the latitude of the Thorodin Dome, British Columbia, *Can. J. Earth Sci.*, **36**(6), 917–943.

Vanderhaeghe, O., Whitney, D., Teyssier, C. & Siddoway, C., 2004. Structural development of the naxos migmatite dome, in *Gneiss Domes in Orogeny*, Vol. **380**, pp. 211–228, eds Whitney, D.L., Teyssier, C. & Siddoway, C.S., Geological Society of America.

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

Vanderhaeghe, O., Guergouz, C., Fabre, C., Duchêne, S. & Baratoux, D., 2019. Secular cooling and crystallization of partially molten Archaean continental crust over 1 Ga, *Comp. Rend. Geosci.*, **351**(8), 562–573.

Vetere, F., Behrens, H., Holtz, F., Vilardo, G. & Ventura, G., 2010. Viscosity of crystal-bearing melts and its implication for magma ascent, *J. Mineral. Petrol. Sci.*, **105**(3), 151–163.

Vidal, M., Gumiaux, C., Cagnard, F., Pouclet, A., Ouattara, G. & Pichon, M., 2009. Evolution of a paleoproterozoic 'weak type' orogeny in the west African Craton (ivory coast), *Tectonophysics*, **477**(3–4), 145–159.

Vigneresse, J.-L., 1995. Crustal regime of deformation and ascent of granitic magma, *Tectonophysics*, **249**(3), 187–202.

Weinberg, F.R. & Schmeling, H., 1992. Polydiapirs: multiwavelength gravity structures, *J. Struct. Geol.*, **14**(4), 425–436.

Weinberg, R.F., 1992. Internal circulation in a buoyant two-fluid newtonian sphere: implications for composed magmatic diapirs, *Earth planet. Sci. Lett.*, **110**(1), 77–94.

Weinberg, R.F., Veveakis, E. & Regenauer-Lieb, K., 2015. Compaction-driven melt segregation in migmatites, *Geology*, **43**(6), 471–474.

Wernicke, B., 1985. Uniform-sense normal simple shear of the continental lithosphere, *Can. J. Earth Sci.*, **22**(1), 108–125.

Whitney, D.L., Teyssier, C. & Vanderhaeghe, O., 2004. Gneiss domes and crustal flow, in *Gneiss Domes in Orogeny*, Vol. **380**, GSA Special Papers, eds Whitney, D.L., Teyssier, C. & Siddoway, C.S., Geological Society of America.

Whittington, A.G., Hofmeister, A.M. & Nabelek, P.I., 2009. Temperature-dependent thermal diffusivity of the Earth's crust and implications for magmatism, *Nature*, **458**(7236), 319–321.

Willet, S.D., 1999. Orogeny and orography: the effects of erosion on the structure of mountain belts, *J. geophys. Res.*, **104**(B12), 28 957–28 981.

Yin, A., 2004. Gneiss domes and gneiss dome systems, *Geol. Soc. Am. Spec. Papers*, **380**, 1–14.

Zibra, I., Lu, Y., Clos, F., Weinberg, R., Peternell, M., Wingate, M., Prause, M., Schiller, M. & Tilhac, R., 2020. Regional-scale polydiapirism predating the Neoarchean Yilgarn Orogeny, *Tectonophysics*, **779**, doi:10.1016/j.tecto.2020.228375.

## APPENDIX A: ADDITIONAL INFORMATION ON THE NUMERICAL ASSUMPTIONS

First, let us recall that the VOF method is well known to preserve material properties at material phases interfaces. Note that in our numerical setup here, while the inclusions can change their volume due to break-up or coalescence, the total mass of each phase is conserved in the domain. In Louis-Napoléon *et al.* (2020, fig. G.1),

we performed tests of mass conservation on multilayer Rayleigh–Taylor problems (including those in Weinberg 1992) and Rayleigh–Benard problems: the mass of each phase was conserved within $10^{-6}$ per cent.

### A1 Complementary description of the Rayleigh numbers for the crustal flow regimes

In Louis-Napoléon *et al.* (2022) we found that distinct flow regimes could be identified depending on two Rayleigh numbers for the unmolten and partially molten domains, $Ra_{UM}$ and $Ra_{PM}$:

$$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 \tilde{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} \tilde{K}_{\text{eff}}^{\text{PM}}}. \quad (A1)$$

In detail, $\Delta T_{UM}$ is the temperature difference between $T_{0\phi}$ at the Earth surface (0 °C) and $T_{Ms}$ at the melting front (here set to 720 °C), $\Delta T_{PM}$ is the temperature difference between $T_{1000\phi}$ at the bottom of the computational domain (1000 °C) and $T_{Ms}$ at the melting front. As for the thickness of the thermal convection regions, and since they are not known *a priori*, we simply take half the depth separating the Earth surface and the bottom of the computational domain, namely $\frac{H_T}{2}$, $\kappa_{UM}$ and $\kappa_{PM}$ are the thermal diffusivities computed with the characteristic mean temperature in each layer, namely $\kappa_{UM} = \kappa((T_{Ms} + T_{0\phi})/2)$ and $\kappa_{PM} = \kappa((T_{Ms} + T_{1000\phi})/2)$, respectively. Finally, we estimate the viscosity using the rheologies described in Section 3.3:

$$\left| \begin{array}{l} \tilde{K}_{\text{eff}}^{\text{UM}} = 0.25 \times 10^6 \times (0.75A)^{-\frac{1}{n}} \times \exp\left(\frac{Q}{nRT_{\text{UM}}}\right), \\ \tilde{K}_{\text{eff}}^{\text{PM}} = \mu_1^0 \times \exp\left[\left(2.5 + \left(\frac{1-M_S}{M_S}\right)^{0.48}\right)(1-M_S)\right], \end{array} \right. \quad (A2)$$

where $\bar{T}_{UM}$ is the characteristic mean temperature in the unmolten layer, $\bar{T}_{UM} = (T_{Ms} + T_{0\phi})/2 = (720 + 0)/2 = 360^{\circ}\text{C}$. In addition, we have $H_T = 45$ km, $\Delta T_{UM} = 720$ K, $\Delta T_{PM} = 280$ K, and $\tilde{K}_{\text{eff}}^{\text{UM}}$ and $\tilde{K}_{\text{eff}}^{\text{PM}}$ depend on the input parameters (cf. eqs 6, 9 and Table 1).

### A2 Simplified diffusivity, heat capacity and latent heat for the temperature equation

Here, $\kappa(T)$ the thermal diffusivity (in m$^2$ s$^{-1}$), and $C_p(T)$ the heat capacity (in m$^2$ s$^{-2}$ K$^{-1}$) are prescribed to depend on temperature according to Whittington *et al.* (2009), in a decoupled manner and similar to, for example Gerya *et al.* (2008):

$$\kappa(T) = \begin{cases} \left(\frac{567.3}{T} - 0.062\right) \times 10^{-6} & \text{if } T < 846\,\text{K}, \\ (0.732 - 1.35 \times 10^{-4}T) \times 10^{-6} & \text{if } T > 846\,\text{K}, \end{cases} \quad (A3)$$

$$C_p(T) = \begin{cases} 902.7 + 0.387T - 22.6 \times 10^{+6}T^{-2} & \text{if } T < 846 \, \text{K}, \\ 1037.6 + 0.146T - 216.7 \times 10^{-6}T^{-2} & \text{if } T > 846 \, \text{K}. \end{cases} \quad (A4)$$

In addition, latent heating is accounted for via the dimensionless apparent heat capacity $C^A$ in the temperature equation. It is defined as $C^A = 1 + \frac{L}{C_p} \frac{\partial M}{\partial T}$, where L is the latent heat (in m$^2$ s$^{-2}$), $C_p$ is the heat capacity (in m$^2$ s$^{-2}$ K$^{-1}$), M is the melt fraction and T is temperature (in K) (Ulvrová *et al.* 2012). Here, we define M using eq. (7), hence we have $\frac{\partial M}{\partial T} = \frac{1}{T_{\text{liq}} - T_{\text{sol}}}$ when $T_{\text{sol}} \le T \le T_{\text{liq}}$ and 0 otherwise. This leads to

$$C^A(T) = 1 + \begin{cases} 0 & \text{if } T < T_{\text{sol}}, \\ St & \text{if } T_{\text{sol}} \le T \le T_{\text{liq}}, \\ 0 & \text{if } T > T_{\text{liq}}, \end{cases} \quad (A5)$$
between horizontal and vertical processes in the formation of the Pilbara Craton, Australia, *Terra Nova*, **19**(1), 1–38.

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

Vanderhaeghe, O. & Teyssier, C., 2001. Partial melting and flow of orogens, *Tectonophysics*, **342**(3–4), 451–472.

Vanderhaeghe, O., Teyssier, C. & Wysoczanski, R., 1999. Structural and geochronological constraints on the role of partial melting during the formation of the Shuswap metamorphic core complex at the latitude of the Thorodin Dome, British Columbia, *Can. J. Earth Sci.*, **36**(6), 917–943.

Vanderhaeghe, O., Whitney, D., Teyssier, C. & Siddoway, C., 2004. Structural development of the naxos migmatite dome, in *Gneiss Domes in Orogeny*, Vol. **380**, pp. 211–228, eds Whitney, D.L., Teyssier, C. & Siddoway, C.S., Geological Society of America.

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

Vanderhaeghe, O., Guergouz, C., Fabre, C., Duchêne, S. & Baratoux, D., 2019. Secular cooling and crystallization of partially molten Archaean continental crust over 1 Ga, *Comp. Rend. Geosci.*, **351**(8), 562–573.

Vetere, F., Behrens, H., Holtz, F., Vilardo, G. & Ventura, G., 2010. Viscosity of crystal-bearing melts and its implication for magma ascent, *J. Mineral. Petrol. Sci.*, **105**(3), 151–163.

Vidal, M., Gumiaux, C., Cagnard, F., Pouclet, A., Ouattara, G. & Pichon, M., 2009. Evolution of a paleoproterozoic 'weak type' orogeny in the west African Craton (ivory coast), *Tectonophysics*, **477**(3–4), 145–159.

Vigneresse, J.-L., 1995. Crustal regime of deformation and ascent of granitic magma, *Tectonophysics*, **249**(3), 187–202.

Weinberg, F.R. & Schmeling, H., 1992. Polydiapirs: multiwavelength gravity structures, *J. Struct. Geol.*, **14**(4), 425–436.

Weinberg, R.F., 1992. Internal circulation in a buoyant two-fluid newtonian sphere: implications for composed magmatic diapirs, *Earth planet. Sci. Lett.*, **110**(1), 77–94.

Weinberg, R.F., Veveakis, E. & Regenauer-Lieb, K., 2015. Compaction-driven melt segregation in migmatites, *Geology*, **43**(6), 471–474.

Wernicke, B., 1985. Uniform-sense normal simple shear of the continental lithosphere, *Can. J. Earth Sci.*, **22**(1), 108–125.

Whitney, D.L., Teyssier, C. & Vanderhaeghe, O., 2004. Gneiss domes and crustal flow, in *Gneiss Domes in Orogeny*, Vol. **380**, GSA Special Papers, eds Whitney, D.L., Teyssier, C. & Siddoway, C.S., Geological Society of America.

Whittington, A.G., Hofmeister, A.M. & Nabelek, P.I., 2009. Temperature-dependent thermal diffusivity of the Earth's crust and implications for magmatism, *Nature*, **458**(7236), 319–321.

Willet, S.D., 1999. Orogeny and orography: the effects of erosion on the structure of mountain belts, *J. geophys. Res.*, **104**(B12), 28 957–28 981.

Yin, A., 2004. Gneiss domes and gneiss dome systems, *Geol. Soc. Am. Spec. Papers*, **380**, 1–14.

Zibra, I., Lu, Y., Clos, F., Weinberg, R., Peternell, M., Wingate, M., Prause, M., Schiller, M. & Tilhac, R., 2020. Regional-scale polydiapirism predating the Neoarchean Yilgarn Orogeny, *Tectonophysics*, **779**, doi:10.1016/j.tecto.2020.228375.

## APPENDIX A: ADDITIONAL INFORMATION ON THE NUMERICAL ASSUMPTIONS

First, let us recall that the VOF method is well known to preserve material properties at material phases interfaces. Note that in our numerical setup here, while the inclusions can change their volume due to break-up or coalescence, the total mass of each phase is conserved in the domain. In Louis-Napoléon et al. (2020, fig. G.1),

we performed tests of mass conservation on multilayer Rayleigh–Taylor problems (including those in Weinberg 1992) and Rayleigh–Benard problems: the mass of each phase was conserved within $10^{-6}$ per cent.

### A1 Complementary description of the Rayleigh numbers for the crustal flow regimes

In Louis-Napoléon et al. (2022) we found that distinct flow regimes could be identified depending on two Rayleigh numbers for the unmolten and partially molten domains, $Ra_{UM}$ and $Ra_{PM}$:

$$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 \tilde{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} \tilde{K}_{\text{eff}}^{\text{PM}}}. \quad (A1)$$

In detail, $\Delta T_{UM}$ is the temperature difference between $T_{0\phi}$ at the Earth surface (0 °C) and $T_{Ms}$ at the melting front (here set to 720 °C), $\Delta T_{PM}$ is the temperature difference between $T_{1000\phi}$ at the bottom of the computational domain (1000 °C) and $T_{Ms}$ at the melting front. As for the thickness of the thermal convection regions, and since they are not known a priori, we simply take half the depth separating the Earth surface and the bottom of the computational domain, namely $\frac{H_T}{2}$. $\kappa_{UM}$ and $\kappa_{PM}$ are the thermal diffusivities computed with the characteristic mean temperature in each layer, namely $\kappa_{UM} = \kappa((T_{Ms} + T_{0\phi})/2)$ and $\kappa_{PM} = \kappa((T_{Ms} + T_{1000\phi})/2)$, respectively. Finally, we estimate the viscosity using the rheologies described in Section 3.3:

$$\begin{cases} \tilde{K}_{\text{eff}}^{\text{UM}} = 0.25 \times 10^6 \times (0.75A)^{-\frac{1}{n}} \times \exp\left(\frac{Q}{nRT_{\text{UM}}}\right), \\ \tilde{K}_{\text{eff}}^{\text{PM}} = \mu_1^0 \times \exp\left[\left(2.5 + \left(\frac{1-M_S}{M_S}\right)^{0.48}\right)(1-M_S)\right], \end{cases} \quad (A2)$$

where $\bar{T}_{UM}$ is the characteristic mean temperature in the unmolten layer, $\bar{T}_{UM} = (T_{Ms} + T_{0\phi})/2 = (720 + 0)/2 = 360^{\circ}\text{C}$. In addition, we have $H_T = 45$ km, $\Delta T_{UM} = 720$ K, $\Delta T_{PM} = 280$ K, and $\tilde{K}_{\text{eff}}^{\text{UM}}$ and $\tilde{K}_{\text{eff}}^{\text{PM}}$ depend on the input parameters (cf. eqs 6, 9 and Table 1).

### A2 Simplified diffusivity, heat capacity and latent heat for the temperature equation

Here, $\kappa(T)$ the thermal diffusivity (in $m^2 s^{-1}$), and $C_p(T)$ the heat capacity (in $m^2 s^{-2} K^{-1}$) are prescribed to depend on temperature according to Whittington et al. (2009), in a decoupled manner and similar to, for example Gerya et al. (2008):

$$\kappa(T) = \begin{cases} \left(\frac{567.3}{T} - 0.062\right) \times 10^{-6} & \text{if } T < 846\text{K}, \\ (0.732 - 1.35 \times 10^{-4}T) \times 10^{-6} & \text{if } T > 846\text{K}, \end{cases} \quad (A3)$$

$$C_p(T) = \begin{cases} 902.7 + 0.387T - 22.6 \times 10^{+6}T^{-2} & \text{if } T < 846 \text{ K}, \\ 1037.6 + 0.146T - 216.7 \times 10^{-6}T^{-2} & \text{if } T > 846 \text{ K}. \end{cases} \quad (A4)$$

In addition, latent heating is accounted for via the dimensionless apparent heat capacity $C^A$ in the temperature equation. It is defined as $C^A = 1 + \frac{L}{C_p} \frac{\partial M}{\partial T}$, where L is the latent heat (in $m^2 s^{-2}$), $C_p$ is the heat capacity (in $m^2 s^{-2} K^{-1}$), M is the melt fraction and T is temperature (in K) (Ulvrová et al. 2012). Here, we define M using eq. (7), hence we have $\frac{\partial M}{\partial T} = \frac{1}{T_{\text{liq}} - T_{\text{sol}}}$ when $T_{\text{sol}} \le T \le T_{\text{liq}}$ and 0 otherwise. This leads to

$$C^A(T) = 1 + \begin{cases} 0 & \text{if } T < T_{\text{sol}}, \\ St & \text{if } T_{\text{sol}} \le T \le T_{\text{liq}}, \\ 0 & \text{if } T > T_{\text{liq}}, \end{cases} \quad (A5)$$
<figure><img src="image_10.png" /><figcaption>**Figure B2.** Time evolution of (a) the melting front depth and (b) the mean vertical velocity, for a series of tests with various initial radii *r* of the inclusions.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure B1.** Zoom on the inclusions in convective cells. Their diameter is 3–4 grid cells, for example about 300–400 m.</figcaption></figure>

where $St = \frac{\mathcal{L}}{C_p(T) \times (T_{\text{liq}} - T_{\text{sol}})}$ is the Stefan number. In this work, $St \sim 1$ since $\mathcal{L} \approx 4 \times 10^5 \text{ m}^2 \text{ s}^{-2}$, $C_p \approx 1200 \text{ m}^2 \text{ s}^{-2} \text{ K}^{-1}$ and $\Delta T = T_{\text{liq}} - T_{\text{sol}} = 400 \text{ K}$.

## APPENDIX B: INFLUENCE OF THE INCLUSIONS SIZE ON THE ‘LAYERING’ REGIME

Here, we test the influence of the inclusions radius *r* on the development of a layer of light inclusions at the top of the convection cell. Therefore, we choose a model case set in the layering regime, hence we reduce further the viscosity pre-factor for the unmolten layer with respect to configuration 1 (cf. eq. 6, Table 2). Namely here, $A_l = 7.81 \times 10^{-6} \text{ Pa}^{-2.3} \text{ s}^{-1}$ (corresponding to case $H_7$ in paper part I). The inclusions radius is varied from *r* = 150 m to *r* = 1200 m. The volume fraction of the inclusions is kept constant, at $\phi \approx 0.35$. The results are presented around time *t* = 10 Myr, which appears as the typical time for dome formation.

In general, we observe that all inclusions entrained in the convective cells progressively adopt a similar equivalent diameter (around 3–4 mesh cells, that is 300–400 m, see Fig. B1). This shows that the size of the inclusions is controlled by the local shear rate, independently of their initial radius.

The temporal evolution of the melting front depth and of the mean vertical velocity are displayed in Fig. B2: it is roughly similar

for all cases. As for the mean vertical velocity, a slow increase occurs during the early stage $t \le 5$ Myr until thermal convection is triggered. After 5 Myr, it increases (together with convection-induced oscillations) and a plateau is reached at about 10 Myr, set in the range [0.25–0.75] cm yr$^{-1}$. We see that for $r \ge 200$ m, the mean vertical velocity becomes roughly independent of *r*. In contrast, for $r \le 150$ m, differences are observed, especially at the onset of convection where the mean velocity is two to three times larger than in the other cases.

Fig. B3 shows that the melting front (red isotherm) reaches a depth of about –25 km in all cases, after 10 Myr. Under the melting front, heterogeneities segregate both upwards (for the light white ones) and downwards (for the black heavy ones), but some also still get dragged into the convection cells. With increasing radius, the segregation of the white layer at ca. 25 km depths becomes more pronounced, while less inclusions remain in the convective cells. The thickness of the white material layer becomes roughly similar (2–3 km) over time, except for $r \le 200$ m, where the formation of the white inclusions layer is not systematic: after checking that mass is conserved, we conclude that the inclusions are undersized with respect to the mesh resolution (a cell size is about 100 m here).

## APPENDIX C: INFLUENCE OF THE INCLUSIONS CONCENTRATION ON THE ‘LAYERING’ REGIME

The concentration of inclusions characterizes the heterogeneity of the partially melting orogenic crust, and we present here tests that highlight its influence on dome formation. First, we display a detailed parametric analysis for one typical model case of the layering regime, and then one specific test for configuration 2 presented in the main text.

### C1 Variation of the inclusions concentration in the ‘layering’ regime case: 2-D simulations

Here, we vary the inclusions concentration ($\phi$, or volume fraction) by changing the distance between inclusions while keeping the inclusions radius fixed, $r = 600$ m. Recall that in 2-D, the theoretical volume fraction of all inclusions (the light and heavy ones)
<figure><img src="image_12.png" /><figcaption>**Figure B2.** Time evolution of (a) the melting front depth and (b) the mean vertical velocity, for a series of tests with various initial radii *r* of the inclusions.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure B1.** Zoom on the inclusions in convective cells. Their diameter is 3–4 grid cells, for example about 300–400 m.</figcaption></figure>

where $St = \frac{\mathcal{L}}{C_p(T) \times (T_{\text{liq}} - T_{\text{sol}})}$ is the Stefan number. In this work, $St \sim 1$ since $\mathcal{L} \approx 4 \times 10^5 \text{ m}^2 \text{ s}^{-2}$, $C_p \approx 1200 \text{ m}^2 \text{ s}^{-2} \text{ K}^{-1}$ and $\Delta T = T_{\text{liq}} - T_{\text{sol}} = 400 \text{ K}$.

## APPENDIX B: INFLUENCE OF THE INCLUSIONS SIZE ON THE ‘LAYERING’ REGIME

Here, we test the influence of the inclusions radius *r* on the development of a layer of light inclusions at the top of the convection cell. Therefore, we choose a model case set in the layering regime, hence we reduce further the viscosity pre-factor for the unmolten layer with respect to configuration 1 (cf. eq. 6, Table 2). Namely here, $A_l = 7.81 \times 10^{-6} \text{ Pa}^{-2.3} \text{ s}^{-1}$ (corresponding to case $H_7$ in paper part I). The inclusions radius is varied from *r* = 150 m to *r* = 1200 m. The volume fraction of the inclusions is kept constant, at $\phi \approx 0.35$. The results are presented around time *t* = 10 Myr, which appears as the typical time for dome formation.

In general, we observe that all inclusions entrained in the convective cells progressively adopt a similar equivalent diameter (around 3–4 mesh cells, that is 300–400 m, see Fig. B1). This shows that the size of the inclusions is controlled by the local shear rate, independently of their initial radius.

The temporal evolution of the melting front depth and of the mean vertical velocity are displayed in Fig. B2: it is roughly similar

for all cases. As for the mean vertical velocity, a slow increase occurs during the early stage $t \le 5$ Myr until thermal convection is triggered. After 5 Myr, it increases (together with convection-induced oscillations) and a plateau is reached at about 10 Myr, set in the range [0.25–0.75] cm yr$^{-1}$. We see that for $r \ge 200$ m, the mean vertical velocity becomes roughly independent of *r*. In contrast, for $r \le 150$ m, differences are observed, especially at the onset of convection where the mean velocity is two to three times larger than in the other cases.

Fig. B3 shows that the melting front (red isotherm) reaches a depth of about –25 km in all cases, after 10 Myr. Under the melting front, heterogeneities segregate both upwards (for the light white ones) and downwards (for the black heavy ones), but some also still get dragged into the convection cells. With increasing radius, the segregation of the white layer at ca. 25 km depths becomes more pronounced, while less inclusions remain in the convective cells. The thickness of the white material layer becomes roughly similar (2–3 km) over time, except for $r \le 200$ m, where the formation of the white inclusions layer is not systematic: after checking that mass is conserved, we conclude that the inclusions are undersized with respect to the mesh resolution (a cell size is about 100 m here).

## APPENDIX C: INFLUENCE OF THE INCLUSIONS CONCENTRATION ON THE ‘LAYERING’ REGIME

The concentration of inclusions characterizes the heterogeneity of the partially melting orogenic crust, and we present here tests that highlight its influence on dome formation. First, we display a detailed parametric analysis for one typical model case of the layering regime, and then one specific test for configuration 2 presented in the main text.

### C1 Variation of the inclusions concentration in the ‘layering’ regime case: 2-D simulations

Here, we vary the inclusions concentration ($\phi$, or volume fraction) by changing the distance between inclusions while keeping the inclusions radius fixed, $r = 600$ m. Recall that in 2-D, the theoretical volume fraction of all inclusions (the light and heavy ones)
<figure><img src="image_12.png" /><figcaption>**Figure B3.** Influence of the inclusions radius on the preservation of the light inclusions layer at ca. 10 Myr (case *H*₁*M*₁*H*₇, fixed *φ* ≈ 0.35). Snapshots display temperature isocontours in Celsius degrees and, on the left (right) half of the domain the viscosity distribution (light and heavy inclusions in white and black, respectively).</figcaption></figure>

is $\phi = \pi r^2/L_i^2$. $\phi$ varies here in the range [0.1; 0.35]. The other model parameters are those of case *H*₇ (cf. Fig. 3), which differs from configuration 1 by the unmolten layer pre-factor viscosity $A_i$ = 7.81 × 10⁻⁶ Pa⁻²s⁻¹ (main text Table 2).

The results are displayed around time 10 Myr in Figs C1 and C2. As in the previous appendix, the inclusions are seen to achieve a roughly similar size due to break-up and coalescence within the convective flow, and independently of the initial volume fraction.
<figure><img src="image_14.png" /><figcaption>**Figure B3.** Influence of the inclusions radius on the preservation of the light inclusions layer at ca. 10 Myr (case *H*₁*M*₁*H*₇, fixed *φ* ≈ 0.35). Snapshots display temperature isocontours in Celsius degrees and, on the left (right) half of the domain the viscosity distribution (light and heavy inclusions in white and black, respectively).</figcaption></figure>

is $\phi = \pi r^2/L_i^2$. $\phi$ varies here in the range [0.1; 0.35]. The other model parameters are those of case *H*₇ (cf. Fig. 3), which differs from configuration 1 by the unmolten layer pre-factor viscosity $A_i = 7.81 \times 10^{-6}$ Pa⁻²s⁻¹ (main text Table 2).

The results are displayed around time 10 Myr in Figs C1 and C2. As in the previous appendix, the inclusions are seen to achieve a roughly similar size due to break-up and coalescence within the convective flow, and independently of the initial volume fraction.
<figure><img src="image_13.png" /><figcaption>**Figure C1.** Time evolution of (a) the melting front depth and (b) the mean vertical velocity, for models in which $\phi$ varies in the range [0.1,0.35]. Here, the inclusions radius $r = 600$ m. Stacking of the light inclusions occurs when concentration exceeds the critical value $\phi = 0.275$.</figcaption></figure>

However, we observe a strong influence of the volume fraction on the development of a light inclusions layer at ca. 20 km depths. Below $\phi = 0.275$, no layer of light inclusions is maintained. This is correlated with a significant change in the average velocity $|U_y|$ between cases with $\phi > 0.275$ and cases with $\phi < 0.275$. This velocity is two to three times larger when $\phi < 0.275$ than when $\phi > 0.275$ (Fig. C1). A transition occurs: below $\phi = 0.275$, convection is barely affected by the inclusions and remains vigorous. The light inclusions are entrained in the convective region. Above $\phi = 0.275$, convection is reduced and a layer of light inclusions can form.

An explanation for the $\phi$-transition stems from considering the thinning of the effective convection layer. Indeed at large $\phi$, the settling of the heavy inclusions observed in Fig. C2 occurs over a greater thickness, and the melting front observed in Fig. C1 is located deeper: both effects act to reduce the thickness of the convective layer. This in turn decreases the Rayleigh number of the unmelted region, reducing the convection intensity and favouring the layering of the light inclusions.

Let us roughly estimate this decrease in Rayleigh number between cases with $\phi = 0.1$ and 0.35. Therefore, we assume that both light and heavy inclusions are stacked into a single layer of length 50 km (that of the computational domain) and of height $h_l^\phi$. The thickness of the convection layer $h_{\text{conv}}^\phi$ can be estimated as $h_{\text{conv}}^\phi = H - h_l^\phi$, where $H$ is the height of the computational domain. In 2-D, we obtain

a ratio of inclusions layer thicknesses $\frac{h_l^{\phi=0.35}}{h_l^{\phi=0.1}} = 3.5$. This leads to a ratio of the convective layer thicknesses $\frac{h_{\text{conv}}^{\phi=0.35}}{h_{\text{conv}}^{\phi=0.1}} = \frac{H-h_l^{\phi=0.35}}{H-h_l^{\phi=0.1}} \approx 0.72$, and thus a ratio between Rayleigh numbers $\frac{Ra^{\phi=0.35}}{Ra^{\phi=0.1}} = (0.72)^3 \approx 0.38$. Hence at high inclusions volume fraction, the local Rayleigh number reduces by more than a factor three and significantly reduces the convection intensity.

Another interpretation might also stem from expressing the effective viscosity of the medium as a Krieger-Dougherty exponential function of $\phi$ (Krieger & Dougherty 1959); however this demonstration requires additional tests that we leave to future studies.

## C2 The influence of $\phi$ on configuration 2 – 3-D simulations

This ultimate test is based on the 3-D simulation *Config$_2$*F$^7-12$ (3-D), where the inclusions concentration is lowered to $\phi = 0.16$. Results Fig C3 are in agreement with observations made in 2-D for low $\phi$ (Figs C1 and C2). At this $\phi$ value below the 0.275 transition, no dome is observed anymore and the rms vertical velocity is two to three times larger than in the reference case (compare with Fig. 9 for $\phi = 0.35$). This 3-D simulation confirms the existence of a $\phi$-transition for domes formation.
<figure><img src="image_15.png" /><figcaption>**Figure C1.** Time evolution of (a) the melting front depth and (b) the mean vertical velocity, for models in which $\phi$ varies in the range [0.1,0.35]. Here, the inclusions radius $r = 600$ m. Stacking of the light inclusions occurs when concentration exceeds the critical value $\phi = 0.275$.</figcaption></figure>

However, we observe a strong influence of the volume fraction on the development of a light inclusions layer at ca. 20 km depths. Below $\phi = 0.275$, no layer of light inclusions is maintained. This is correlated with a significant change in the average velocity $|U_y|$ between cases with $\phi > 0.275$ and cases with $\phi < 0.275$. This velocity is two to three times larger when $\phi < 0.275$ than when $\phi > 0.275$ (Fig. C1). A transition occurs: below $\phi = 0.275$, convection is barely affected by the inclusions and remains vigorous. The light inclusions are entrained in the convective region. Above $\phi = 0.275$, convection is reduced and a layer of light inclusions can form.

An explanation for the $\phi$-transition stems from considering the thinning of the effective convection layer. Indeed at large $\phi$, the settling of the heavy inclusions observed in Fig. C2 occurs over a greater thickness, and the melting front observed in Fig. C1 is located deeper: both effects act to reduce the thickness of the convective layer. This in turn decreases the Rayleigh number of the unmelted region, reducing the convection intensity and favouring the layering of the light inclusions.

Let us roughly estimate this decrease in Rayleigh number between cases with $\phi = 0.1$ and 0.35. Therefore, we assume that both light and heavy inclusions are stacked into a single layer of length 50 km (that of the computational domain) and of height $h_l^\phi$. The thickness of the convection layer $h_{\text{conv}}^\phi$ can be estimated as $h_{\text{conv}}^\phi = H - h_l^\phi$, where $H$ is the height of the computational domain. In 2-D, we obtain

a ratio of inclusions layer thicknesses $\frac{h_l^{\phi=0.35}}{h_l^{\phi=0.1}} = 3.5$. This leads to a ratio of the convective layer thicknesses $\frac{h_{\text{conv}}^{\phi=0.35}}{h_{\text{conv}}^{\phi=0.1}} = \frac{H-h_l^{\phi=0.35}}{H-h_l^{\phi=0.1}} \approx 0.72$, and thus a ratio between Rayleigh numbers $\frac{Ra^{\phi=0.35}}{Ra^{\phi=0.1}} = (0.72)^3 \approx 0.38$. Hence at high inclusions volume fraction, the local Rayleigh number reduces by more than a factor three and significantly reduces the convection intensity.

Another interpretation might also stem from expressing the effective viscosity of the medium as a Krieger-Dougherty exponential function of $\phi$ (Krieger & Dougherty 1959); however this demonstration requires additional tests that we leave to future studies.

## C2 The influence of $\phi$ on configuration 2 – 3-D simulations

This ultimate test is based on the 3-D simulation *Config$_2$*F$^7-12$ (3-D), where the inclusions concentration is lowered to $\phi = 0.16$. Results Fig C3 are in agreement with observations made in 2-D for low $\phi$ (Figs C1 and C2). At this $\phi$ value below the 0.275 transition, no dome is observed anymore and the rms vertical velocity is two to three times larger than in the reference case (compare with Fig. 9 for $\phi = 0.35$). This 3-D simulation confirms the existence of a $\phi$-transition for domes formation.
<figure><img src="image_14.png" /><figcaption>**Figure C2.** Influence of the inclusions concentration $\phi$ ($H_1M_1H_7$, $r = 600$ m). Snapshots at $ca. 10$ and $15$ Myr: for $\phi > 0.275$, the white inclusions layer remains despite convection underneath. Temperature isocontours in Celsius degrees and, on the left (right) half of the domain the viscosity distribution (light and heavy inclusions in white and black, respectively).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure C2.** Influence of the inclusions concentration $\phi$ ($H_1M_1H_7$, $r = 600$ m). Snapshots at $ca. 10$ and $15$ Myr: for $\phi > 0.275$, the white inclusions layer remains despite convection underneath. Temperature isocontours in Celsius degrees and, on the left (right) half of the domain the viscosity distribution (light and heavy inclusions in white and black, respectively).</figcaption></figure>
<figure><img src="image_15.png" /></figure>

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

<figure><img src="image_17.png" /><figcaption>**Figure C3.** Configuration 2 with a lower concentration of inclusions, $\phi = 0.16$ (and here, $t_1 = 10$ Myr and $t_2 = 15$ Myr). (a) Time evolution of the flow with 3-D views and 2-D vertical profiles along plane $z = 25$ km (see legend of Fig. 5). (b) 3-D trajectory of some markers with initial location $(x, y, z) = (25, -40, 25)$, $(35, -22, 25)$ and $(40, -37, 30)$ marked by circles. (c) Markers depth location, temperature and velocity over time.</figcaption></figure>
<figure><img src="image_17.png" /></figure>

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

<figure><img src="image_19.png" /><figcaption>**Figure C3.** Configuration 2 with a lower concentration of inclusions, $\phi = 0.16$ (and here, $t_1 = 10$ Myr and $t_2 = 15$ Myr). (a) Time evolution of the flow with 3-D views and 2-D vertical profiles along plane $z = 25$ km (see legend of Fig. 5). (b) 3-D trajectory of some markers with initial location $(x, y, z) = (25, -40, 25)$, $(35, -22, 25)$ and $(40, -37, 30)$ marked by circles. (c) Markers depth location, temperature and velocity over time.</figcaption></figure>