# Porous Melt Flow in Continental Crust-A Numerical Modeling Study

P. Maierová, P. Hasalová, Karel Schulmann, P. Štípská, O. Souček

► To cite this version:

P. Maierová, P. Hasalová, Karel Schulmann, P. Štípská, O. Souček. Porous Melt Flow in Continental Crust-A Numerical Modeling Study. Journal of Geophysical Research: Solid Earth, 2023, 128 (8), 10.1029/2023JB026523. [https://doi.org/10.1029/2023JB026523](https://doi.org/10.1029/2023JB026523)

HAL Id: hal-04221346

https://hal.science/hal-04221346v1

Submitted on 29 Sep 2023

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2023JB026523

**Key Points:**

* Numerical thermo-mechanical models simulate porous melt flow through hot felsic continental crust

* Melt can propagate through pores with the spacing of 1 mm or larger with a velocity of millimeters to centimeters per year and form large partially molten zones

* Style of melt flow and solid deformation vary between melt-enhanced convection, diapirism, and porosity waves depending on model parameters

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

P. Maierová,
petra.maierova@geology.cz

**Citation:**

Maierová, P., Hasalová, P., Schulmann, K., Štípská, P., & Souček, O. (2023). Porous melt flow in continental crust—A numerical modeling study. *Journal of Geophysical Research: Solid Earth*, 128, e2023JB026523. https://doi.org/10.1029/2023JB026523

Received 6 FEB 2023

Accepted 19 JUL 2023

**Author Contributions:**

**Conceptualization:** P. Maierová, P.

Hasalová, K. Schulmann, P. Štípská

**Methodology:** P. Maierová, O. Souček

**Software:** P. Maierová, O. Souček

**Supervision:** P. Hasalová, K. Schulmann,

P. Štípská, O. Souček

**Validation:** P. Štípská, O. Souček

**Visualization:** P. Maierová

**Writing—original draft:** P. Maierová, P. Hasalová, K. Schulmann, P. Štípská,

O. Souček

© 2023. The Authors.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

## Porous Melt Flow in Continental Crust—A Numerical Modeling Study

P. Maierová¹ ID, P. Hasalová¹, K. Schulmann¹,², P. Štípská¹ ID, and O. Souček³

¹Centre for Lithospheric Research, Czech Geological Survey, Prague, Czech Republic, ²EOST, Institute de Physique du Globe, Université de Strasbourg, Strasbourg, France, ³Faculty of Mathematics and Physics, Mathematical Institute, Charles University, Prague, Czech Republic

**Abstract** In continental crust, rapid melt flow through macroscopic conduits is usually envisaged as the most efficient form of melt transport. In contrast, there is growing evidence that in hot continental crust, grain-scale to meso-scale porous melt flow may operate over long distances and over millions of years. Here, we investigate the dynamics of such porous melt flow by means of two-dimensional thermo-mechanical numerical models using the code ASPECT. Our models are crustal-scale and describe the network of pores through which the melt flows by permeability that depends on the spacing of the pores. Our results suggest that assuming realistic material properties, melt can slowly migrate in the hot and thick continental crust through pores with a characteristic spacing of 1 mm or larger. Despite its low velocity (millimeters to centimeters per year), over millions of years, such flow can create large partially molten zones in the middle-lower crust and significantly affect its thermal state, deformation, and composition. We examined the role of the permeability, melt and solid viscosities, the slope of the melting curve and temperature conditions. We obtained contrasting styles of melt distribution, melt flow, and solid deformation, which can be categorized as melt-enhanced convection, growth of partially molten diapirs and melt percolation in porosity waves. Our numerical experiments further indicate that grain-scale porous flow is more likely in rocks where the melt productivity increases slowly with temperature, such as in metaigneous rocks.

**Plain Language Summary** Continental crust contains a significant proportion of rocks that can melt relatively easily. Geophysical methods that measure, for example, seismic velocities or electric conductivity revealed large partially molten regions in the continental crust. In addition, solidified melt pathways (veins, dykes) are often observed in rocks that were formerly in hot continental interiors. Usually, melt flow through macroscopic pathways is assumed to be the most efficient type of melt transport. However, several geological studies have shown that melt may flow through much finer (microscopic) pores, which form a network along boundaries of rock grains. The question arises, how efficient this microscopic melt flow is, and how it affects the properties of the continents. We address this issue using numerical models. According to our models, microscopic melt flow is very slow—melt moves millimeters to centimeters per year. Despite its low velocity, over millions of years, it forms large partially molten zones, significantly increases the temperature in the crust, changes its composition and softens it. The style of this process varies notably depending on the rock properties and temperature conditions. For example, in some models, kilometer-sized melt batches move through nearly immobile rock, while in others, melt and rock both move together as a stirred mush.

## 1. Introduction

The common view of melt migration in the continental crust involves melt segregation from the original sites of melt generation by grain-scale melt flow and then flow through a semi-to interconnected network of leucosomes (Figure 1a), veins/dykes or ductile fractures (e.g., Petford et al., 2000; Sawyer, 2001; Weinberg & Regenauer-Lieb, 2010) into ascent conduits that ultimately feed plutons at shallow crustal levels (e.g., Brown, 2010; Sawyer, 2001). Networks marked by leucosomes and/or veins and dykes were described in both metasedimentary and metaigneous migmatitic terranes worldwide (e.g., Hall & Kisters, 2012; Korhonen et al., 2010; Tanner, 1999; Yakymchuk et al., 2013). Melt flow through leucosomes is envisaged as the dominant mechanism of mesoscale melt migration (Brown, 2010; Vanderhaeghe, 1999).

In contrast to this view, there is growing geological evidence for an entirely different style of melt migration in hot continental crust. Several examples were found in the eastern European Variscan belt in the Bohemian Massif,
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Two contrasting types of migmatites resulting from different styles of melt migration: (a) migmatite with foliation-parallel leucosomes representing former well-defined channels of melt migration. Example from Zanskar, NW Himalaya. (b–c) Migmatite without any obvious melt segregation channels as a result of porous melt flow. Example from Eger crystalline unit, NW Bohemian Massif. (c) Sketch of melt distribution in migmatite without obvious melt segregation.</figcaption></figure>

where the denuded middle and lower crust of an ancient hot orogen crop out (e.g., Schulmann et al., 2009). A peculiar package of mid-crustal metaigneous gneisses and migmatites (~10,000 km²) was described (Hasalová, Janoušek, Schulmann, Štípská, & Erban, 2008; Hasalová et al., 2008a, 2008b) where leucosomes are rare and former melt is distributed in the whole rock volume along grain boundaries (Figures 1b and 1c). The origin of these migmatites cannot be explained by in situ melting only and requires pervasive melt migration from an external, although chemically similar, source. The overall amount of passing melt was large, but the melt fraction at any time was small, significantly lower than required for the breakdown of the solid crystal framework (~25%–40% of melt, e.g., Rosenberg & Handy, 2005). Besides these mid-crustal migmatites, a similar process was identified several hundred kilometers away in exhumed felsic lower crust, which underwent subduction (Aguilar et al., 2020; Štípská et al., 2019; Závada et al., 2018, 2021). In addition, U-Pb monazite geochronology shows that the interaction between the melt and the solid rock was a long-lasting process that could have taken several millions of years (8–15 Myr; Závada et al., 2021). All this suggests that in these migmatites, melt segregation was inefficient and grain-scale pervasive melt flow was an important mechanism of melt transport throughout the crust. This process was described by Hasalová et al. (2008a) as a grain-scale, diffuse and reactive flow of a felsic silicate melt through a rock, when melt passes pervasively along grain boundaries through the whole rock volume and changes the macro- and microscopic appearance of the rock, as well as its mineral and whole-rock composition. A similar process is well documented in mantle rocks (e.g., Kelemen et al., 1992), where flow channels are recognized by contrasting mineralogy and composition to their surroundings (e.g., dunite in peridotite). Adopting the terminology commonly used for mantle rocks, we refer to this process as porous melt flow also for crustal rocks.

There are indications that melt migrated in this manner also in other places. In the Vosges Mountains, central European Variscan Belt in France, a large portion of partially molten crust (granitoids, orthogneisses/migmatites) was infiltrated pervasively by granitic melt that changed their macro- and microscopic appearance and composition (Hasalová et al., 2015). In Fiordland, New Zealand, a large arc region (>10,000 km²) of lower-crustal gabbroic gneissic rocks was interpreted to record the flow of granitic melt through microscopic pathways along grain boundaries (Daczko et al., 2016; Meek et al., 2019; Stuart et al., 2016, 2018). Here, melt migrated through large-scale shear zones that formed channels where hydrous melts reacted with the host rock. In these rocks, reaction microstructures that resulted from porous melt flow were described. In the Karakoram Shear Zone, NW Himalaya, an extensive dyke swarm, where dykes form an interconnected network was described (Reichardt & Weinberg, 2012). These dykes acted as magma pathways to feed large magma bodies, but each dyke records internal porous melt flow (Hasalová et al., 2011).

A number of key questions arise from these observations: under what conditions can porous melt flow operate in the continental crust? How does micro-scale porous flow affect the properties of the crust in comparison with the
meso-macro-scale flow in conduits? And what is the actual style of melt migration and the resulting melt distribution in the interiors of hot orogens? These questions need to be answered as the melt transport and distribution have profound implications on the evolution of the continental crust: (a) the low viscosity of partially molten rocks affects the internal deformation of orogens and their topography (e.g., Jamieson et al., 2011; Vanderhaeghe & Teyssier, 2001); (b) heat transported by melt affects temperature distribution (e.g., Leitch & Weinberg, 2002) and (again) viscosity and deformation of crustal rocks (e.g., Brown, 2010; Hyndman, 2019; Sandiford & McLaren, 2002); and (c) partial melting and melt migration controls chemical differentiation of the crust and stabilization of continents (e.g., Hawkesworth & Kemp, 2006; Rudnick et al., 2003; Sawyer et al., 2011). The melt distribution in active orogens has been partly depicted by geophysical studies that have shown the presence of magmatic zones (MZs) and partially molten regions in the crust (e.g., overview by Unsworth, 2010). However, the resolution of geophysical methods is not sufficient to reveal the melt distribution at a smaller scale, and discrimination between the contrasting mechanisms of melt transport is therefore impossible (e.g., Pritchard & Gregg, 2016).

Further insight into melt-migration processes in the continental crust can be provided by numerical modeling. Thermal effects of melt transport were studied using models where the melt flow and/or site of emplacement were prescribed kinematically (e.g., Annen et al., 2015; Leitch & Weinberg, 2002). Besides that, some thermo-mechanical models of crustal evolution contain kinematically prescribed melt flow (Cao et al., 2016) or approximations of the processes of melt segregation and emplacement (e.g., Babeyko et al., 2002). However, dynamics of melt migration in the continental crust has been modeled only rarely. Most of the published numerical models of melt migration focus on the mantle and therefore do not take into account specific properties of the crust (e.g., Connolly & Podladchikov, 2007; Dohmen & Schmeling, 2021) or neglect thermal evolution (Keller et al., 2013). The notable exception is one-dimensional (1D) models of porous flow coupled with thermal evolution and melting that were used to investigate melt migration in the lower crust induced by basaltic intrusions (e.g., Jackson et al., 2005; Solano et al., 2014). In these models, it is assumed that the partially molten rock consists of a solid porous matrix and melt that fills the pores entirely. The melt fraction is therefore identical with the porosity of the rock. These models show the formation of buoyancy-driven porosity waves (PWs, regions of higher porosity that move with respect to the solid) (Spiegelman, 1993) and moving high-porosity layers. In these regions, the melt is transported over several kilometers in total in less than a few millions of years. Importantly, Solano et al. (2014) were able to simulate chemical variations between the solid rock and passing melt due to compositional equilibration between them. The limiting factor of these 1D models is that they cannot describe the lateral component of, for example, buoyancy-driven motion.

Recent models by Schmeling et al. (2023, 2019) overcome this shortcoming as they simulate porous flow, thermal evolution, and melting in two dimensions (2D). The models by Schmeling et al. (2019) were mainly designed to explain the formation of the magmatic body in the middle crust beneath the Altiplano plateau (Ward et al., 2014) as a result of high heat flux at the bottom of the crust and high initial temperature. They reproduced an ascent of a large partially molten diapir from the hot lower crust leading to the formation of a high-porosity “magmatic” body in the middle crust over a time span of tens to hundreds of thousands of years. The main focus of their study was on the porous flow in the upper “magmatic” part of the diapir with high melt fractions, while the evolution of the regions with less than ~20% of melt was discussed only marginally. The very recent study by Schmeling et al. (2023) further explored a similar model setup, but assuming several heat pulses at the bottom of the crust. They obtained two main modes of model evolution, which they called “batholith emplacement” and “convective recycling” modes. While the batholith emplacement mode is characterized by the formation of high-porosity “magmatic” bodies, lower porosities (<20%) are typical for the convective recycling mode. In summary, dynamic numerical models have confirmed the viability of porous melt flow in the crust. However, only a few models investigated the range of conditions typical for the migmatites that witness micro-scale porous melt flow: mid-to high-temperature, mid-to lower-crustal pressure, felsic host rock, granitic melt composition, and most importantly, melt fraction below the threshold required for disaggregation of the solid crystal framework.

Among the possible reasons for the lack of relevant numerical studies are simplified theoretical considerations suggesting that the porous melt flow is inefficient in the continental crust (e.g., Rutter et al., 2006). Their main argument is the high viscosity of silicate melt (several orders of magnitude higher than the viscosity of the basaltic melt [e.g., Giordano & Dingwell, 2003]) causing a less efficient flow of the silicate melt through the solid rock matrix (McKenzie, 1984). While this general conclusion is valid, experiments show that the silicate hydrous melt viscosity can be as low as ~10⁴ Pa s under mid-crustal conditions and possibly even lower in the lower crust
(Holtz et al., 2001; Weinberg & Hasalová, 2015; Whittington et al., 2004). With such a viscosity, grain-scale melt migration cannot be a-priori neglected, particularly if the crust is maintained at supra-solidus conditions for millions of years (e.g., Cavalcante et al., 2018; Wolfram et al., 2019).

Here, we incorporate the current knowledge of crustal melting and properties of crustal rocks and granitic melt into 2D numerical models of melt migration in hot orogens. We evaluate the efficiency of the porous melt flow for different material properties and thermal conditions within the bounds given by experiments and geological and geophysical observations. We use a similar approach as Schmeling et al. (2019), but we assume long-lasting high-temperature conditions at the base of the crust. This mimics the high thermal regime typical in some active orogens as a result of, for example, thinned mantle lithosphere (Hyndman, 2019). First, we briefly describe the model setup including the reasoning behind the choice of individual model parameters (Section 2). Then the modeling results are described in Section 3. Finally, we discuss the results in terms of possible styles of melt migration and development of partially molten zones, and we compare them with observations from hot orogens, as well as with the previously published numerical models (Section 4). In the end, we conclude that numerical modeling identifies porous melt flow as a viable, efficient, and long-lasting mechanism of melt transport in continental crust, and suggests that it may represent an important melt transfer process.

## 2. Numerical Model Setup

In this section, we give a brief overview of the model setup. More details on the form of the solved equations and the choice of material parameters are in Text S1 in Supporting Information S1. The solved equations were implemented in the ASPECT code version 2.3.0 (Bangerth, Dannberg, Gassmoeller, & Heister, 2020; Bangerth, Dannberg, Gassmöller, et al., 2020; Dannberg & Heister, 2016; Heister et al., 2017; Kronbichler et al., 2012) that was used for calculation of the presented results. The melting parametrization, treatment of material composition, and definition of some material properties (density, viscosity, and permeability) were modified from the standard definition in ASPECT, as described in Text S1 in Supporting Information S1.

We aim to study porous melt flow as a mechanism of large-scale melt transport in felsic continental crust. Our models are crustal scale and their duration is several millions of years. All the properties that are below the model resolution, which is hundreds of meters in the partially molten regions, are averaged. The continental crust is described as a material consisting of two phases—solid rock matrix and liquid melt. Both phases can deform as well as move relative to each other, and transformation between them corresponds to melting and solidification. The volumetric fraction of the melt is one of the model variables, and because it is equal to the volume fraction of pores in the matrix, it is also called porosity.

We solve the equations of conservation of mass and momentum of the solid and melt phases that are both considered incompressible. The conservation of momentum in the melt is approximated by the Darcy law:

$$ \mathbf{u}_f - \mathbf{u}_s = -\frac{K_D}{\phi}(\nabla p_f - \rho_f \mathbf{g}), \quad (1) $$

where $\mathbf{u}_s$ and $\mathbf{u}_f$ are velocities of the solid and melt, respectively, $K_D$ is the ratio of the permeability and melt viscosity called the Darcy coefficient, $\phi$ is the porosity, $p_f$ is the melt pressure, $\rho_f$ is the melt density and $\mathbf{g}$ is the gravity. We solve the equation of heat conservation taking into account heat advection by the solid and melt, diffusion of heat, volumetric (radiogenic) heating and latent heat of melting.

Melting is parametrized using the solidus and liquidus temperatures (Figure 2a) that depend linearly on the composition c and lithostatic pressure P:

$$ \begin{aligned} T_S &= T_{S0} + \Delta T_P P + \Delta T_c c, \\ T_L &= T_S + \Delta T_{SL}. \end{aligned} \quad (2) $$

The composition c is calculated from the compositions of the solid and melt: $c = c_s(1 - \phi) + c_f\phi$. These compositions are advected by the velocities of the solid and melt, respectively, and they are adjusted during melting/solidification so that they lie on the solidus and liquidus lines, respectively (Figure 2a). From the melting parametrization, it follows that the composition c can be regarded as a measure of depletion in constituents that preferentially enter the melt during partial melting. In the description of the model results, we use the terms
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Melting parametrization. In equilibrium, the compositions of the solid and liquid phases lie on the solidus and liquidus lines for a given pressure, temperature, and bulk composition. The distance between the solidus and liquidus lines is given by the $\Delta T_{SL}$ parameter. (b) Sketch of nonlinear melting curves for two contrasting rock types with a high and low water content (rock A and B, respectively). The melting curve can be approximated by a linear function for a given temperature. The slope of the linear approximation then defines the $\Delta T_{SL}$ parameter.</figcaption></figure>

“depleted” for the material composition of $c > c_0$ and “enriched” for the composition of $c < c_0$, where $c_0 = 0$ is the reference composition. For interpretation of the model results, the composition c can be regarded as an approximate measure of the silica content and of the “dryness” of the rock (the silica and water contents decrease with increasing c). Model materials with strongly depleted and enriched compositions are proxies for the felsic granulite and granite, respectively.

The adopted melting parametrization is highly simplified and it is valid only in a limited range of conditions and compositions. For felsic crustal rocks, the reference solidus $T_{50}$ and the pressure dependence of solidus $\Delta T_p$ are relatively well constrained and we fix their values to 650°C and 100 K GPa$^{-1}$, respectively. The solidus-liquidus difference $\Delta T_{SL}$ is related to the local slope of the melting curve (Figure 2b). This property differs significantly depending on the rock type, and we test different values of $\Delta T_{SL}$ in the range of 200°C–1000°C (e.g., Clemens, 2006; Clemens & Vielzeuf, 1987).

We model warm and weak lower and middle crust where we assume that viscous creep dominates deformation, and we neglect brittle deformation. We are aware that brittle deformation may occur even under these conditions, especially because the pressure of melt reduces the effective pressure in the rock. We quantify this possibility using the compaction pressure $p_c$, which is the same as the effective pressure in partially molten material (Keller et al., 2013). A negative compaction pressure is related to melt overpressure and is associated with decompaction and matrix dilatation. In places where the melt overpressure overcomes the rock cohesion, tensile brittle failure may occur. In our considerations, we take the value of 5 MPa as an estimate of the rock cohesion and show where $p_c < -5$ MPa in our models. In these places, tensile failure might occur, followed by the formation of dykes. These processes are, however, out of the scope of this study.

We assume that the density depends on temperature and composition with the reference densities of the solid and melt of 2,800 and 2,400 kg m$^{-3}$, respectively (e.g., Christensen & Mooney, 1995; Whittington et al., 2004). For the shear viscosity of the solid matrix, we use a modified power law:

$$ \eta_s = \eta_0 \dot{\epsilon}_{II}^{\frac{1-n}{n}} \exp\left(\frac{A}{nRT}\right) \exp(\alpha_c c_s) \exp(-\alpha_{\phi} \phi), \quad (3) $$

where $\dot{\epsilon}_{II}$ is the strain rate, the parameters $\eta_0$, $n$, and $A$ are given by the quartzite flow law (Ranalli, 1995) scaled 1×, 0.5×, and 0.25× in different model cases (Figure 3a, solid, dashed, and dotted black lines), the parameter $\alpha_c$ determines the dependence on the composition and the parameter $\alpha_{\phi}$ determines the porosity-dependent weakening (cf. Rudge, 2018). The viscosity predicted by the adopted flow law is ~10$^{19}$ Pa s for the hot unmolten crust and it is further reduced in the presence of melt to ~10$^{18}$–10$^{16}$ Pa s for $\phi = 0.1–0.3$ (Figure 3b, solid black line) (cf. Clark & Royden, 2000; Dong et al., 2020). The value of 10$^{16}$ Pa s is also the lower bound imposed on the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Shear and bulk viscosities of the solid matrix. (a) Selected experimental and theoretical flow laws for the shear viscosity of felsic crustal rocks for $\dot{\epsilon}_{II} = 10^{-14} s^{-1}$ and along the initial temperature profile used in our models. Red—granite (Ranalli, 1995); green—quartzite (Lu & Jiang, 2019; Ranalli, 1995, and references therein); pink—felsic granulite (Wen et al., 2021). Dry quartzite flow law by Ranalli (1995) (solid black line) is adopted in this study. Weaker flow laws (dashed and dotted gray) are obtained by reducing the viscosity 2× and 4×. (b) Dependence of shear and bulk viscosity on the porosity, $\eta_s \sim \exp(-\alpha_\phi\phi)$ and $\eta_b \sim \exp(-\alpha_\phi\phi)(\beta_1 + \beta_2 \frac{1}{\phi})$ for $\alpha_\phi = 20$, $\beta_1 = 2$ and $\beta_2 = 0.3$ (solid and dashed black lines). Theoretical dependence after Schmeling et al. (2012) for elliptic melt inclusions with the aspect ratio of 0.1 and 0.03 (solid and dashed dark and light gray lines) is shown for comparison.</figcaption></figure>

shear viscosity in order to prevent rapid small-scale convection and related numerical issues. Even though the solid material itself is incompressible, the porous solid matrix is compressible and has a finite bulk viscosity. We define the bulk viscosity as the shear viscosity multiplied by a factor that blows up like $1/\phi$ for $\phi \to 0$ (Figure 3b, dashed black line) (e.g., Schmeling et al., 2012; Simpson et al., 2010). The viscosity of the melt is assumed to be constant, $\eta_f = 10^4$ Pa s (cf. Holtz et al., 2001), and its role in the model evolution can be tested by changing the Darcy coefficient (see below).

The permeability $k_\phi$ depends on the porosity, the spacing of the pores $d$ and their geometry. In our models, we adopt the definition of permeability proposed by Keller et al. (2013):

$$k_{\phi} = k_{0}d^{2}\phi^{3}(1 - \phi)^{2}. \qquad (4)$$

The prefactor $k_0$ is related to the geometry of the pores, and its value has been estimated mostly between $10^{-1}$ and $10^{-3}$ (e.g., Connolly et al., 2009; Wark & Watson, 1998). The spacing of the pores $d$ is usually assumed to be equal to the grain size, that is $\sim 1$ mm, because the melt resides along the grain boundaries. However, in a rock with a system of dykes or other melt-rich zones, $d$ may rather correspond to the spacing of these zones (Schmeling et al., 2019), and with different values of $d$, we may approximate migration of melt through such a system. In the solved equations, the permeability and the melt viscosity occur only in the Darcy coefficient. We, therefore, define the reference Darcy coefficient:

$$K_{D0} = \frac{k_0 d^2}{\eta_f}. \qquad (5)$$

We test values between $10^{-12}$ and $10^{-9} m^2 Pa^{-1} s^{-1}$, where the intermediate value of $10^{-11} m^2 Pa^{-1} s^{-1}$ corresponds to the grain size of 1 mm, melt viscosity of $10^4$ Pa s and $k_0$ of $10^{-1}$, which means that it is probably a maximum value for grain-scale porous flow in crustal rocks.

The model domain (Figure 4) represents a cross-section through thickened felsic continental crust: the top boundary corresponds to the surface of the

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Model domain and boundary conditions (BC). The model domain is 50 km wide and 50 km deep. Gray lines show an example of the adaptively refined computational mesh (model REF at $t = 2$ Myr). At the bottom boundary, the temperature is increased by a cosine anomaly schematically depicted by the blue line.</figcaption></figure>
Earth and the bottom boundary is located in the lower crust. We assume that below this level, there is no crustal material that could melt. The domain is rectangular, 50 km wide and 50 km deep. The assumed large crustal thickness is based on data from collisional orogens (e.g., Beck & Zandt, 2002; Schulmann et al., 2009; Zhang et al., 2011). Boundary conditions for the velocity are: free slip at the bottom boundary, zero velocity at the sides of the domain and zero traction (i.e., mechanically open) top boundary. The porosity is initially zero everywhere. The initial composition is zero except in the lowermost 1 km of the domain, where we prescribe a depleted zone that precludes melting next to the boundary and possible grid-scale numerical oscillations of the porosity. The initial temperature is given by an analytical solution of heat conduction for T = 0°C at the top boundary and depth-dependent radiogenic heat sources. Near the bottom boundary, the initial temperature is slightly above 800°C, which is just below the solidus temperature in that region. During the time evolution, zero heat flux is prescribed on the sides, T = 0°C at the top boundary and the initial temperature with an added high-temperature anomaly at the bottom boundary (Figure 4). The temperature anomaly lasts during the first 5 Myr of model evolution and then it is switched off. The amplitude of the anomaly $T_a$ is one of the model parameters.

The finite element method is used for spatial discretization and solution of the equations (Dannberg & Heister, 2016). The finite elements are quadrilateral, quadratic in the velocities, temperature, porosity and compositions, and linear in the pressures. The mesh resolution is variable (adaptively refined), with the smallest and largest possible element size of ~390 m and ~6 km, respectively.

# 3. Results

## 3.1. Reference Model

First, we describe the evolution of the reference model. The reference model has intermediate values of material parameters: The reference Darcy coefficient is set to $10^{-11} \text{m}^2 \text{Pa}^{-1} \text{s}^{-1}$. The shear viscosity of the solid is prescribed by the quartzite flow law scaled by 0.5 (dashed black line in Figure 3). Both the solidus-liquidus difference and the amplitude of the temperature anomaly are prescribed to be 300°C. This particular reference model setup was selected because it represents a case where melt migrates efficiently along the grain boundaries, and because it shows the most important features that are, to varying degrees, also present in other model cases: partially molten and magmatic zones, porosity waves, significant solid deformation and large- and small-scale compositional variations. The parameters of the reference model are summarized in Table S1 in Supporting Information S1. Melt fractions, isotherms, and velocities in representative time steps are shown in Figure 5. The entire time evolution of the reference model is available in Movie 1.

Initially, the temperature in the entire model domain is below the solidus and there is no melt present. As the temperature anomaly at the bottom boundary supplies heat, temperature at the bottom of the domain quickly increases and melting occurs. Gradually, a diapir-shaped region with ~10%-20% of melt develops ("partially molten zone," PMZ, red to yellow color in Figure 5, left column). In the initial stage, a cap with more than 30% melt ("magmatic zone," MZ) is formed on top of the PMZ (white color in Figure 5, left column, 0.5 Myr). The coexistence of PMZ and MZ lasts approximately 1 Myr. Afterward, the MZ rises to a shallower level where it cools down and solidifies (Figure 5, 2 Myr). Melt in the PMZ is spatially organized into so-called PWs—circular or planar regions of increased melt fraction (which is the same as porosity in the two-phase solid-melt material) that move with respect to the solid material. In this model, the PWs are kilometer-sized and merge into channels along the sides of the PMZ. This style of melt migration lasts several millions of years, although the average melt fraction in the PMZ decreases (Figure 5, 4 and 6 Myr). Melt percolation continues after the temperature anomaly at the base has been switched off at 5 Myr of model evolution. At 7 Myr of model evolution, that is 2 Myr after the switch off, nearly all the melt has solidified.

The melt moves through the solid at a maximum relative velocity $u_{rel}$ of up to ~5 cm yr$^{-1}$, decreasing with time (Figure 5, right column, colors), with typical values around 1 cm yr$^{-1}$. The magnitude of $u_{rel}$ depends mainly on the porosity due to the porosity dependence of permeability (Equation 4). The motion of the melt through the solid is nearly vertical. In addition to the flow of the melt with respect to the solid, the solid material deforms (Figure 5, right column, white arrows). The deformation in the PMZ is facilitated by weakening of the solid material due to its porosity (Equation 3). Thanks to this weakening, stresses induced by buoyancy forces are sufficient to drive visible solid deformation. A single upwelling—a diapir—centered above the temperature anomaly and downwellings along its sides develop. The diapiric upwelling is also partly responsible for the uprise of the MZ
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Four representative time steps of the reference model (model REF with $K_{DO} = 10^{-11} m^2 Pa^{-1} s^{-1}$, solid viscosity given by the quartzite flow law (Ranalli, 1995) scaled 0.5x, $\Delta T_{SL} = 300°C$ and $T_a = 300°C$). Tics on x- and y-axes are at every 10 km. Left: melt fraction (which is the same as the porosity), isotherms at every 100°C from 500°C to 1100°C, and the -5 MPa isoline of compaction pressure (green contour, present only at 0.5 Myr time step in this model). Right: velocity of solid (white arrows) and magnitude of $u_{rel}$, the relative velocity of melt with respect to the solid (colors in the background). The direction of $u_{rel}$ is always nearly vertically upward.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Composition *c* in the reference model at the end of the simulation when all the melt has solidified. (a) Negative values (pink) and positive values (gray) correspond to enrichment and depletion, respectively, in the component that preferentially enters the melt. (b) Composition profile at *x* = 25 km. The impact of composition on solidus and liquidus temperatures is given by the Δ*T*<sub>c</sub> term (Equation 2). In the strongly enriched middle crust, the decrease of the melting temperature due to enrichment reaches ~120°C. In the lowermost crust, the melting temperature increases by 300°C.</figcaption></figure>

and its cooling. The velocity of the solid is ~2 cm yr⁻¹ at the maximum, and in the PMZ, it is lower or comparable in magnitude to *u*<sub>rel</sub>. In the late stage, without any melt, the solid motion ceases due to a lower temperature at the bottom boundary and higher viscosity even though the diapiric body remains buoyant. Throughout the simulation, the melt overpressure in the PMZ is too low to facilitate brittle failure. The regions with melt overpressure exceeding 5 MPa are limited to the margin of the MZ (Figure 5, 0.5 Myr).

The motion of the solid has a dominant effect on the temperature field because the velocities of the melt and the solid are comparable, and the melt fraction is mostly below 20%. The rise of the hot material causes an increase in temperature at mid-crustal depth (see isotherms in Figure 5, left column). In the diapir at a 30–40 km depth, the temperature increases by ~100°C–150°C compared to the initial temperature. The latent heat released during solidification contributes to this increase mainly in the uppermost part of the PMZ and MZ, where most of the solidification takes place. The temperature increase due to latent heat release can be estimated as φΔST/ *c*<sub>p</sub>, where ΔS is the entropy change due to melting, *T* is temperature and *c*<sub>p</sub> is the heat capacity. As the average melt fraction at 30–40 km depth is less than ~10%, the increase of temperature due to the latent heat release is below ~30°C.

During the model evolution, a compositionally enriched middle crust and depleted lower crust gradually form (Figure 6 shows the distribution of composition at the end of model evolution). As a result of this enrichment-depletion (and of the increase of temperature in the middle crust) the PMZ moves from the lower crust toward the middle crust. The solidus of the enriched material defines the top of the PMZ and the solidus of the depleted material defines its bottom in the lower crust. While the lower crust shows an almost vertical gradient of composition, in the middle crust, there are kilometer-scale compositional variations that develop from PWs. These variations are initially circular, but they become flattened and distorted due to the vertical shortening of the head of the diapir. Due to the composition-dependent solidus temperature, the enriched remnants of PWs stay partially molten for a longer time than their surroundings and form “fingers” at the top margin of the PMZ (Figure 5, left column, 4 and 6 Myr).

The compositional variations affect the density and viscosity of the material. In the lower crust, the density increases due to depletion by up to 100 kg m⁻³ and the shear viscosity increases by ~1 order of magnitude. In the enriched middle crust, a complementary decrease in density and viscosity is observed. Anyway, the effect of composition (*c*<sub>s</sub>) on the viscosity is only minor compared to the porosity-induced weakening. The solid viscosity in the PMZ decreases due to the porosity by 1–2 orders of magnitude to the value of ~10¹⁸ Pa s. In the centers of PWs, it is lowered to ~10¹⁷ Pa s and in the MZ, it drops to 10¹⁶ Pa s, which is the lowermost cut-off value imposed in our models.

## 3.2. Parametric Study

We will now describe the effect of different parameters on model evolution using a series of models with varied parameter values. First, we focus on the effect of the Darcy coefficient and the solid viscosity. Then, we focus on the role of melting parameters and temperature conditions. Model parameters are the same as in the reference model except for those directly specified in the text and figures. The full list of model cases is in Table S2 in Supporting Information S1. Representative time steps for selected models are presented in Figures 7 and 8. Typical styles of model evolution are summarized in Figure 9.

### 3.2.1. Darcy Coefficient and Viscosity of Solid

Darcy coefficient *K*<sub>D</sub>, which is the ratio of permeability and melt viscosity, is the key parameter governing the relative motion between the melt and the solid matrix. In models with a high Darcy coefficient (model KD-10 with *K*<sub>D0</sub> = 10<sup>-10</sup>m² Pa⁻¹ s⁻¹, first row in Figure 7) melt percolates more easily and the relative velocity *u*<sub>rel</sub> is ~4× higher than in the reference model. Consequently, the typical porosity and thus melt-related buoyancy is lower and, as a result, the solid velocity is ~2× smaller than in the reference model. Slower solid motion results in a smaller temperature perturbation (max. ~100°C). The melt is clearly focused in PWs that gather into vertical
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Results for four selected models: high Darcy coefficient (model KD-10 with $K_{D0} = 10^{-10} m^2 Pa^{-1} s^{-1}$), low Darcy coefficient (model KD-12 with $K_{D0} = 10^{-12} m^2 Pa^{-1} s^{-1}$), low viscosity (model V025 with solid viscosity given by the quartzite flow law [Ranalli, 1995] scaled 0.25×) and high viscosity (model V1 with solid viscosity given by the quartzite flow law [Ranalli, 1995]). Tics on x- and y-axes are at every 10 km. Left: Melt fraction (which is the same as the porosity), isotherms at every 100°C from 500°C to 1100°C, and the -5 MPa isoline of compaction pressure (green contour). Middle: velocity of solid (white arrows) and magnitude of $u_{rel}$, the relative velocity of the melt with respect to the solid (colors in the background). The direction of $u_{rel}$ is always nearly vertically upward. Right: composition at the end of the simulation. Negative values (pink) and positive values (gray) correspond to enrichment and depletion, respectively, in the component that preferentially enters the melt.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Results for four selected models: low $\Delta T_{SL}$ and high $T_a/\Delta T_{SL}$ ratio (model DTSL200Ta220 with $\Delta T_{SL} = 200^{\circ}$C, $T_a = 220^{\circ}$C, that is $1.1 \times \Delta T_{SL}$), high $\Delta T_{SL}$ and high $T_a/\Delta T_{SL}$ ratio (model DTSL500Ta550 with $\Delta T_{SL} = 500^{\circ}$C, $T_a = 550^{\circ}$C, that is $1.1 \times \Delta T_{SL}$), low $\Delta T_{SL}$ and low $T_a/\Delta T_{SL}$ ratio (model DTSL200Ta100 with $\Delta T_{SL} = 200^{\circ}$C, $T_a = 100^{\circ}$C, that is $0.5 \times \Delta T_{SL}$), high $\Delta T_{SL}$ and low $T_a/\Delta T_{SL}$ ratio (model DTSL1000Ta500 with $\Delta T_{SL} = 1000^{\circ}$C, $T_a = 500^{\circ}$C, that is $0.5 \times \Delta T_{SL}$). The figure layout is the same as in Figure 7.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Summary of melt migration styles in the different modeled cases during the first 5 Myr of evolution. C—melt-enhanced convection, D—growth of a partially molten diapir, PW—melt percolation in porosity waves. Numbers in italics denote the time when melt (or solidified melt) first rises to the 30 km depth (Myr). Yellow dot—models with magmatic zones with more than ~30%–40% of melt (small = until 0.5 Myr, mid = until 2 Myr, big = until 5 Myr). Green dot—models with compaction pressure locally less than –5 MPa (small = until 0.5 Myr, mid = until 2 Myr, big = until 5 Myr). The reference model is in boldface. (a) Models with different reference Darcy coefficient $K_{D0}$ and solid viscosity $\eta_s$ (scaled quartzite flow law; Ranalli, 1995). (b) Models with different $\Delta T_{SL}$ and $T_a/\Delta T_{SL}$ ratio. Gray contour lines show absolute values of temperature anomalies.</figcaption></figure>

channels or “fingers.” The migrating melt alters the composition in these channels toward enrichment which lowers the solidus temperature, which enhances the focusing of melt into the channels. The focusing mechanism is then a combination of mechanically driven melt flow directed into the channel and compositionally enhanced melt production in the channel. The tips of the partially molten fingers rise to a shallower depth than the MZ in the reference model, and as a result of the composition-induced lowering of the solidus, they reach lower temperatures. The rapid transport of melt into the relatively cold middle crust causes solidification without the formation of any MZ, and an overall lower melt fraction in the PMZ. The solidified melt forms several separated lobes of enriched material. The magnitude of enrichment in these lobes (and complementary depletion in the lower crust) is higher than in the reference model. The strong enrichment causes unrealistically low solidus temperature in the tips of the lobes and the presence of melt at $T < 600^{\circ}$C. The fast melt percolation into the cold and high-viscosity regions induces a higher melt overpressure at the top margin of the PMZ (green contours in Figure 7, left panel). Inside the PMZ, the magnitude of the compaction pressure remains small.

In the opposite case with a low Darcy coefficient (model KD-12 with $K_{D0} = 10^{-12} m^2 Pa^{-1} s^{-1}$, the second row in Figure 7), melt flow through the solid matrix is hindered. As a result, the melt is retained near the site of melting in the lower crust. As the melt fraction in this region increases, the solid matrix becomes weaker and vigorous melt-enhanced convection occurs. During the first 0.5 Myr, the melt fraction is around 30% and quickly decreases with time to ~20% as the lowermost crust gradually becomes more depleted. The fluid pressure generated by melting is not sufficient to induce tensional fractures neither in the PMZ nor at its contact with the colder crust above it. This is because the melt percolation into the cold crust and associated decompaction of this crust are very limited. The depleted lower crust and enriched middle crust form, but the maxima of depletion and enrichment are significantly smaller than in the reference model. Small-scale compositional variations develop in the PMZ, but they are efficiently stirred by convection.

The solid shear viscosity $\eta_s$ directly affects the deformation of the solid matrix. In addition, in our formulation, it is linked with the bulk viscosity of the solid, which affects the volumetric compaction and decompaction of the matrix and, therefore, the melt flow. In the models with low solid viscosity (model V025 with the quartzite flow law scaled 0.25x, third row in Figure 7), the solid material easily responds to the stresses induced by buoyancy. That results in efficient convection even if the melt fraction is relatively low (<20%). Several (two) convective cells develop, their shapes change with time and the typical solid velocity is 2–3× higher in this model than in the reference case. This vigorous convection stirs the PMZ and the nearby unmolten material. With the melt fraction below 20% and the melt velocity comparable to or smaller than the solid velocity, most of the heat is advected by the solid. Due to the high solid velocity, the top of the PMZ reaches a shallower depth than in the reference model, and the temperature perturbation in the middle crust is more elevated. The relative velocity $u_{rel}$ is similar
to or slightly smaller than in the reference model. As a result, enriched and depleted materials separate quite efficiently and, despite some mixing, the enriched middle crust and depleted lower crust develop. The magnitude of the compaction pressure is generally low in models with a low solid viscosity.

A high solid viscosity (model V1 with the quartzite flow law, fourth row in Figure 7) prevents the development of convection and hinders the growth of the partially molten diapir. PWs move through an otherwise almost static solid, merge and form a MZ. As there is negligible upwelling of the solid, the PMZ remains at deeper levels than in the reference model and the temperature perturbation in this region is <100°C. The high viscous strength of the material increases the potential for brittle failure. Still, the sites with $p_c < -5$ MPa are located on top of the MZ or near the sides of the PMZ and not within the PMZ. The final compositional anomalies are more pronounced than in the reference model.

### 3.2.2. Temperature Anomaly and Solidus-Liquidus Difference

In the second series of models, we varied the difference between the solidus and liquidus temperatures ($\Delta T_{SL}$ parameter) and the temperature anomaly ($T_a$). The change of the $\Delta T_{SL}$ changes the liquidus temperature, while the solidus for a given pressure and reference composition remains the same (see Equation 2). In models with a larger $\Delta T_{SL}$, the same temperature variation has a smaller effect on the melt fraction compared to the models with a smaller $\Delta T_{SL}$. It is then the ratio of $T_a$ and $\Delta T_{SL}$ that governs the amount of melt in each model. Besides that, the $\Delta T_{SL}$ affects the compositional dependence of the solidus and liquidus temperatures: As we keep the $\Delta c_{SL}$ parameter fixed, the change of $\Delta T_{SL}$ affects the slope of both the solidus and liquidus lines (see Figure 2a). Consequently, a high $\Delta T_{SL}$ implies a strong dependence of the solidus and liquidus temperatures on the composition and vice versa. The enrichment-depletion patterns that develop in the models therefore substantially modify the melting temperatures in the high $\Delta T_{SL}$ models.

In models with a high $T_a/\Delta T_{SL}$ ratio (models DTSL200Ta220, DTSL500Ta550, first and second rows in Figure 8), the temperature at the bottom boundary (initial temperature plus temperature anomaly) is significantly higher than the solidus temperature at the bottom of the model domain for the reference composition. A lot of melt is therefore created in the lower crust during the first time steps and the melt quickly percolates upwards and forms a MZ on top of the PMZ. The MZ is relatively large, but it solidifies quickly due to the diapiric upwelling and its small remnants are present only in model DTSL500Ta550 at 2 Myr. In both these models, the high melt fraction induces the weakening of the material and fast diapiric or convective motion. In the cooler model DTSL200Ta220, the PMZ is narrower and convection develops only shortly at ~2.5–3 Myr. A clear depletion-enrichment pattern is observed in these models. In contrast, models with lower $T_a/\Delta T_{SL}$ ratio show a moderate melt fraction (around 10%) that percolates in PWs and the absence of MZs (models DTSL200Ta100, DTSL1000Ta500, third and fourth rows in Figure 8).

In models with a large $\Delta T_{SL}$ (model DTSL1000Ta500, fourth row in Figure 8), a wide PMZ forms that extend to shallow levels and low temperatures. Its vertical extent results from the strong compositional dependence of the solidus temperature that induces melting of the enriched material in the middle crust despite its relatively low temperature. The upper boundary of the PMZ shows several lobes that have an enriched composition. At the tips of these lobes, sites with a high melt overpressure are located. In contrast, for a small $\Delta T_{SL}$ (e.g., model DTSL200Ta100, third row in Figure 8), the melt is spatially focused in a narrow diapir-shaped zone. In line with that, compositional variations at the end of the model evolution show a round upper margin without lobes (right column in Figure 8).

### 3.2.3. Trends in Model Characteristics

In addition to Figures 5–8, we summarize the results of the full series of models in Figure 9. For each model case, the following characteristics are indicated: the presence of MZs and the possibility of brittle deformation, the time scale of melt uprise and the dominant style of melt migration and solid deformation in the PMZ—melt-enhanced convection (C), growth of partially molten diapirs (D) or melt percolation in PWs. By melt-enhanced convection, we denote a process where the buoyancy of the melt and the hot solid efficiently drives the motion of the solid. As a result, the solid material moves rapidly (several centimeters per year) and forms several upwellings and downwellings whose shapes change with time. The partially molten diapir is also driven by buoyancy, but its shape is rather static and only gradually grows in vertical and lateral directions. PWs are regions of higher melt fraction that move with respect to the solid. They are mostly circular in our models, and they may focus into channels or “fingers” of higher pulsating porosity. For the formation of PWs, movement of melt through the solid is essential.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** A schematic picture of the main observed styles of model evolution: melt-enhanced convection (C), diapirism (D), and melt percolation in porosity waves (PW). Red shading depicts the melt fraction. PWs (spherical “bubbles” of increased melt fraction) transfer also melt inside the partially molten diapir, as sketched in the middle panel. The dotted line is an isotherm showing that the highest temperature increase in the middle crust occurs in models with convection (left panel). On the other hand, melt percolation in PWs (right panel) can bring melt to cooler levels of the middle crust. Striped areas depict possible magmatic zones. They are formed mainly in models with a high-temperature (high $T_a/\Delta T_{SL}$ ratio) and therefore a high melt fraction. The green dots and arrows show sites with a potential for brittle failure as a result of high (>5 MPa) melt overpressure.</figcaption></figure>

These three styles observed in our models are schematically depicted in Figure 10. In some models, several of these styles are present at the same time or develop consecutively. For example, diapirism in the early stage of evolution may change to convection and then, in the late stage, to PWs (e.g., in model DTSL500Ta500). Usually, PWs are visible inside the partially molten diapirs and, to a certain degree, also in the models where convection dominates.

The dominant style of deformation and melt migration (see “C,” “D,” and “PW” in Figure 9) depends on the velocities of the melt and solid that in turn depend mainly on the Darcy coefficient and solid viscosity, respectively. Models with a high $K_D$ and high $\eta_s$ (top right corner in Figure 9a) develop melt percolation in PWs. In models with a low $K_D$ and low $\eta_s$, diapiric or convective motion develops (bottom left in Figure 9a). In the second series of models, the style of deformation varies mainly due to the temperature and porosity dependence of the viscosity. PWs dominate in cold and, therefore, high-viscosity models (models with a small temperature anomaly, bottom left corner in Figure 9b), while diapirs or convection develop in high-temperature models. The solid and melt velocities further determine the characteristic time scale of the model evolution. For each model, we evaluated the time it takes until the first melt (or solidified melt) reaches the depth of 30 km, that is until it rises 20 km above the bottom boundary (numbers in italics in Figure 9). We observe that this time increases with the solid viscosity and decreases with $K_D$ (longest times are in the top left corner in Figure 9a). In the second series of models, this time mainly decreases with the temperature anomaly (in Figure 9b, longest times are in the bottom left corner, where absolute values of $T_a$ are lowest).

Green dots in Figure 9 indicate the potential for brittle failure characterized by the presence of fluid overpressure above 5 MPa in the models. The potential for brittle failure increases with $K_D$, with the solid viscosity (toward top right in Figure 9a), and with $\Delta T_{SL}$ because with a large $\Delta T_{SL}$, melt percolates into relatively cold high-viscosity material (right in Figure 9b). Yellow dots in Figure 9 indicate whether MZs develop in the models showing that they are present in models with a high solid viscosity (top in Figure 9a) and a high-temperature anomaly (top in Figure 9b). The models with MZs and with a high fluid overpressure mostly coincide. In such cases, melting on top of the MZs might potentially induce brittle fracturing and these fractures could then form pathways for the flow of melt toward the surface. These trends are in general agreement with predictions based on the analysis of the solved equations, as elaborated in detail in Text S2 in Supporting Information S1. Let us note, however, that our estimates of the compaction pressure stem from employing purely viscous matrix rheology. Considering more complex rheology involving additional relaxation mechanisms (e.g., visco-elasticity) would probably lead to a
reduction of the compaction pressure in the marginal zones. From this point of view, our interpretations of the potential for brittle failure should be approached with caution.

## 4. Discussion

In this section, we discuss the modeled styles of melt flow and their implications in terms of melt distribution, the associated solid matrix deformation and composition. We compare our results with geophysical and geological data from modern and ancient orogens. We point out important model limitations and compare our results with previous numerical-modeling studies.

### 4.1. Melt Transport Styles and Development of Partially Molten and Magmatic Zones in Models

The presented models suggest that the gravity-driven porous melt flow might be an important mechanism of melt transport in hot continental crust. The melt flow velocity ($u_{rel}$ in our models) depends mainly on the Darcy coefficient, or more specifically on the rock permeability that depends on the spacing of the melt-filled pores. For the spacing of the pores ~1 mm, corresponding to melt distributed at grain boundaries, the melt percolates through the rock at the velocity of millimeters to centimeters per year. For a higher permeability, corresponding to meso-macro-scale spacing of pores, the melt velocities are significantly higher. Conversely, for a lower Darcy coefficient (e.g., due to the spacing of pores smaller than 1 mm), the melt is effectively locked within the rock. Melt-induced weakening then facilitates the deformation of the rock while the melt is carried along. While this basic relationship is clear already from the analysis of the equations of porous flow, our models show its impact on melt distribution and flow on the crustal scale.

In our models, the melt flow is closely related to the development of PMZs that are dominant features of the model lower and middle crust (Figure 10). These zones are large, up to ~25 km thick, typically contain ~5%–25% volume fraction of melt, and their character depends on the model parameters. Large PMZs where melt percolates in PWs develop in models with relatively small temperature anomalies and gradual melting curves (low $T_a/\Delta T_{SL}$ ratio), high Darcy coefficient and high solid viscosity. A typical feature of these zones is a relatively stable melt fraction that oscillates only due to the passage of PWs. In contrast, in models with a low Darcy coefficient, low solid viscosity and high melt fractions (high $T_a/\Delta T_{SL}$ ratio), diapirism or convection dominate. In these models, there is enhanced solid matrix deformation in the PMZs, regions with different melt fractions mix and in some cases, the surrounding unmolten material may be entrained into the PMZ.

On top of the PMZs, MZs (striped regions in Figure 10) form. They form mainly in models with a moderate-high Darcy coefficient and solid viscosity, and a high typical melt fraction (high $T_a/\Delta T_{SL}$ ratio). It means that for the same temperature conditions, MZs form more likely in material with a steep melting curve (small $\Delta T_{SL}$; Figure 2b). In contrast, formation of long-lasting PMZs with the melt fraction around 10% and PWs is preferred in material with more gradual melting curves (large $\Delta T_{SL}$; Figure 2b). The values of the compaction pressure in our models suggest that brittle fracturing is unlikely inside the PMZ or MZ, but may occur at their margins (green dots in Figure 10). Brittle fracturing would therefore most likely occur at the contact of melt with cold and strong crust. In a real crust, brittle fracturing (both tensile and shear) in these regions would lead to the draining of the melt from the top of the molten zones into the overlying sub-solidus crust.

### 4.2. Geophysical Evidence for Partially Molten and Magmatic Zones in Hot Orogens

Our models show the development of large and long-lasting PMZs accompanied by less prominent and transient MZs. The presence of such partially molten and MZs in the continental crust of active hot orogens is evidenced by geophysical studies. In the Himalayan-Tibetan system, seismic “bright spots” (e.g., Gaillard et al., 2004) and low electrical resistivity of ~1–3 Ωm (Wei et al., 2001) suggest that MZs form lenticular bodies that are typically present at a ~15–20 km depth beneath South Tibet. In Central Tibet, the zone of low resistivity deepens and in some regions continues even into the mantle (Wei et al., 2001). The resistivity in this zone shows higher values (~5–10 Ωm) consistent with the presence of ~10% of melt. Many more geophysical studies suggest that MZs in the middle crust and partially molten regions in the middle-lower crust are present in active orogens such as the Himalayas, Anatolian Plateau, or Andes (e.g., overview by Hacker et al., 2014; Unsworth, 2010). Locations of
MZ and PMZ and melt fractions in our models are in broad agreement with these observations, although the MZs form deeper in our models due to the relatively low temperature in the model middle crust.

Recently, a high-resolution magnetotelluric study by Dong et al. (2020) revealed significant spatial variations of conductivity beneath Tibet. They resolved ~100 km spaced vertical zones of high conductivity surrounded by relatively low-conductivity crust. The calculated vertically averaged melt fractions varied from ~15% to 25% in the high-conductivity zones to ~5%–10% of melt in the surrounding region. Based on our models, we speculate that this pattern may result from the diapiric/convective motion of the partially molten crust. Similar structures, interpreted as partially molten diapirs, were imaged beneath the Altiplano plateau using gravity data (del Potro et al., 2013). These structures were ~15 km wide, rooted at a 20 km depth, and their melt fraction varied from zero (fully crystallized granitoid rock) to almost 100% melt in their cores occupying less than 5% of the volume.

All these geophysical observations suggest that large partially molten zones are common in hot continental crust. This view is further supported by geochemical and geological data from migmatitic and magmatic crustal rocks that witness the state of fossil partially molten and magmatic zones (Schwindinger & Weinberg, 2017; Weinberg et al., 2021). An important piece of evidence is long residence times (up to millions of years) of crystals in contact with melt (e.g., Cavalcante et al., 2018) that contrast with short lifetimes of magma chambers (e.g., Bachmann & Bergantz, 2004) and direct field observations related to the interaction between flowing melt and solid framework (for a summary see Weinberg et al., 2021). Altogether, during the last decade, these observations caused a fundamental shift in the understanding of crustal melting that now accentuates the role of partially molten rocks with a solid crystal framework—so-called mushes. It is now proposed that mush zones can occupy not only the lower crust (e.g., Annen et al., 2015; Delph et al., 2017), but that they may form large interconnected transcrustal columns (Cashman et al., 2017; Sparks et al., 2019). The porous melt flow (both micro- and meso-scale) has been considered an important process in crustal mush zones (e.g., Sparks et al., 2019; Weinberg et al., 2021). In this perspective, the PMZs that develop in our models due to the porous melt flow correspond to these transcrustal mush zones. Our models confirm that PMZs can be dominant and long-lasting features in hot continental crust, while MZs develop only in some cases and may be short-lived. In addition, the models show a significant variability in the shapes, deformation and thermal structure of the PMZs.

## 4.3. Geological Evidence for Porous Melt Flow

Here, we will compare our modeling results with the main features of the lower-to-mid-crustal migmatites that once formed a part of ancient partially molten zones in hot orogens and that were affected by porous melt flow. We compare the modeled and observed melt fractions, discuss the relationship between the rock type (protolith) and melt flow style, changes in rock composition due to chemical equilibration, the time scale of porous flow and development of brittle deformation.

### 4.3.1. Melt Fraction

In our models, the instantaneous melt fraction in the PMZ is ~5%–25%. The melt fraction is higher in models where convection or diapirism develop (~15%–20% of melt), while it is smaller in PMZs where the melt is transported mainly by PWs (~10% of melt). In most models, melt fraction is spatially variable in the PMZ—in the centers of the PWs and in channels that form at the sides of the PMZ, the melt fraction is 5%–10% higher than average.

It is tempting to compare these calculated melt fractions with the values of melt fractions inferred from geological data. However, accurate quantification of the amount of melt in migmatites remains difficult. While we can estimate relatively well the amount of former melt crystallized in the rock based on microstructural and field observations (e.g., Hasalová et al., 2011; Sawyer, 2001; Závada et al., 2018), it is impossible to accurately determine the instantaneous fraction of melt that was present in the rock and that corresponds to the melt fraction calculated in our models as well as to the one estimated by geophysical methods.

Reported crystallized-melt fractions in migmatitic terranes showing grain-scale porous flow are 5%–10% (e.g., Hasalová et al., 2008a). At higher fractions, around 15% of melt, leucosomes begin to form (Stuart et al., 2018). The instantaneous melt fraction was constrained by Schulmann et al. (2008) and Závada et al. (2007, 2018), who investigated the deformational characteristics of migmatites with up to 15% of crystallized melt. They concluded that these rocks preserved the solid skeleton and the instantaneous melt fraction therefore has not reached the
threshold required for its disaggregation (e.g., Rosenberg & Handy, 2005). For migmatites with more than ~20% of melt, a well-defined leucosome network is typical (e.g., Vernon & Paterson, 2001; Yakymchuk et al., 2013). In some of these cases, thermodynamic modeling showed that the crystallized-melt fractions significantly exceeded the amount of melt that could be locally produced in the rock. This indicates that external melt migrated into and/ or through these migmatites and crystallized there (Štípská et al., 2019; Yakymchuk et al., 2013).

The instantaneous melt fraction in the PMZs in our models (~5%–25%) is below the disaggregation threshold. In this respect, our models agree well with the observations. The modeled melt fraction is typically somewhat higher than the crystallized-melt fraction in migmatites with a grain-scale network of pores; the higher values are comparable to the melt fractions in leucosome-rich migmatites. Still, some models show porous flow in PMZ with the melt fraction around 10% (e.g., model DTSL1000Ta500), which is similar to the crystallized-melt fraction measured in migmatites showing grain-scale porous flow.

### 4.3.2. Rock Type—Metaigneous Versus Metasedimentary Protoliths

In our models, the rock type is characterized by its melting parameters and viscosity. The melting parametrization is extremely simplified and we focused only on the effect of one parameter, $\Delta T_{SL}$, that is related to the slope of the melting curve (Figure 2). Models with slowly increasing melting curves (large $\Delta T_{SL}$) are characterized by low melt fractions, large PMZs and melt transport in PWs (e.g., Figure 8, model DTSL1000Ta500). Only for very high temperature conditions (high $T_g$), we obtain higher melt fractions, formation of MZs and important solid deformation—diapirism or convection. In models with steep melting curves (small $\Delta T_{SL}$), high melt fractions can be reached at moderately high temperatures (e.g., Figure 8, model DTSL200Ta220), and these models develop diapirism or convection and MZs. The style of evolution is further affected by the solid rock viscosity: in high-viscosity material, the solid deformation is slower and the porous melt flow is relatively more important compared to low-viscosity models. The high-viscosity models also show overall higher melt fractions and larger compositional variations. These trends can be compared with the characteristics of the two distinct lithologies of felsic continental crust—metaigneous and metasedimentary rocks.

Felsic metaigneous rocks are relatively dry and usually contain less than 10% of hydrous minerals. As a result, their melting curve increases slowly with increasing temperature and at crustal high-temperature conditions, they are able to produce only around 10% of melt (~20% in extreme cases) (Hasalová et al., 2008c; Štípská et al., 2019; Yakymchuk et al., 2013). So far, it is only in metaigneous migmatites where the grain-scale porous melt flow has been detected (Štípská et al., 2019, and references therein). In these rocks melt fraction generally remains distributed throughout the rock volume, although in some places, melt has locally accumulated in dilatant or compressive structures such as shear-bands, or foliation (e.g., Závada et al., 2018). However, melt drainage networks are not widely described. Similarly, interstitial melt flow in crystallizing mush in magma chambers has been described (Weinberg et al., 2021). These observations agree with our models showing that for a slowly increasing melting curve (large $\Delta T_{SL}$), which approximates melting in metaigneous rocks, porous melt flow in regions with a relatively small melt fraction (~10%) is typical.

In contrast, metasedimentary rocks are rich in hydrous minerals and in hot lower to middle crust may produce up to ~50% of melt, which is far more than average metaigneous rocks. The range of possible melt fractions strongly depends on the chemical composition of the protolith (e.g., Clemens, 2006; le Breton & Thompson, 1988; Patiño Douce & Harris, 1998; White et al., 2001). The typical appearance of metasedimentary migmatite is stromatite—migmatite with a visible melt-drainage network (Figure 1a). This suggests that melt segregates easily in these rocks and is drained via the segregation network. Distribution of the leucosomes is controlled by deformation and rock anisotropy (e.g., Hall & Kisters, 2012; Vernon & Paterson, 2001; Yakymchuk et al., 2013) that is, together with very low viscosity, typical for these rocks. In natural rocks, segregation of melt into leucosomes may be also enhanced by the fast and in some cases abrupt increase of melt fraction over a limited temperature interval. The melt flow through the grain-boundary network then might be too slow to balance the melt production. Heterogeneity and anisotropy of these rocks then lead to the formation of melt-rich and melt-poor zones that will eventually transform into a meso-scale network of melt conduits (e.g., Brown, 2004).

The abrupt melting and melt segregation that are characteristic for metasediments cannot be well approximated in our models. The models with steep melting curves (small $\Delta T_{SL}$) typically show high melt fractions that induce weakening of the solid matrix and substantial solid deformation. However, as the models do not simulate the melt segregation into meso-scale conduits, we can not directly compare these results with the melting process
in metasedimentary rocks. In an evolved stage of melting in metasedimentary rocks, we can assume that the meso-scale network of pores is already formed throughout the rock, and we can approximate this situation by a model with high permeability (i.e., high Darcy coefficient). The high-permeability models evolve at shorter time scales, the melt fraction in the PMZ is smaller and PWs dominate melt transport (e.g., model KD-10, Figure 7). This also results in lower temperatures in the PMZ compared to the models with lower permeability. MZs are spatially limited to the tops of the lobes of the PMZ. However, we have not tuned the model to contain both the high permeability and the steep melting curve, the two typical features of metasedimentary rocks.

### 4.3.3. Melt-Host Rock Chemical Equilibration

Chemical equilibration between host rock and melt is an important process during porous melt flow (e.g., Hasalová et al., 2008c; Meek et al., 2019). In our models, it is approximated by the equilibration of the composition c. The adopted melting parametrization (Figure 2a), migration of melt through the solid and their compositional equilibration altogether result in the formation of depleted (c > 0) and enriched (c < 0) regions in the model. Melt migration and composition equilibration between melt and solid lead to the differentiation between depleted model lower crust and enriched model middle crust. The enriched model middle crust shows compositional variations at a kilometer-scale that develop from regions of a different melt fraction, such as the PWs (e.g., Figure 6).

In continental crust, the newly formed melt is typically enriched in incompatible and mobile elements (e.g., Rudnick et al., 2003; Sawyer et al., 2011). This melt then starts to migrate and along its migration path, it partially crystallizes and/or equilibrates with the surrounding rocks. As a result, the host rock is enriched in elements carried by the passing melt. Upward movement of the melt then drives differentiation between depleted lower crust and granitic upper crust (e.g., Rudnick et al., 2003). Chemical heterogeneities are observed in outcrop- to kilometer-scale where they reflect different degrees of equilibration between the rock and the passing melt (Meek et al., 2019; Stuart et al., 2018).

The large-scale crustal differentiation and formation of compositional heterogeneities are to a certain degree mimicked by our models, which in addition suggest that the kilometer-scale chemical variations in migmatitic terranes may result from a variable amount of melt that migrated through the rock. The outcrop-scale variations are below the model resolution. In real rocks, these variations cannot result from PWs (they have different characteristic sizes), and they probably develop from inherited heterogeneities enhanced by interaction with the migrating melt.

### 4.3.4. Timescale of “Migmatization”

The models show that with a sufficient and persistent heat source (a prescribed temperature anomaly in our case), the melt flow in the PMZ may be active for millions of years. If the heat source does not change, the typical melt fraction in the PMZ gradually decays due to progressive differentiation of the crust leading to depletion of the hot lower crust. Once the heat source is switched off, the PMZ can survive for up to one to 2 millions of years (Figure 5). The duration of the PMZ after the switch off of the heat source depends on the actual condition imposed at the bottom boundary. In our model setup, we instantly lower the temperature at the bottom boundary, which implies an unrealistic heat outflow from the domain. A more natural boundary condition resulting in a gradual decrease in the temperature at the bottom boundary would prolong the duration of the PMZ in the middle crust.

The time scales of our models agree with the longevity of partially molten zones in hot orogens. For example, migmatites from different parts of highly metamorphosed units in the Himalayas suggest that partial melting took up to 25 Myr (e.g., Ding et al., 2021). The ancient Araçuaí orogen in Brazil remained partially molten for at least 25 Myr (Cavalcante et al., 2018). The Famatinian orogeny in NW Argentina shows a 60 Myr record of cyclic melting (e.g., Wolfram et al., 2019). Long duration of high-temperature conditions (hundreds of thousands of years to millions of years) was proposed for mush zones (e.g., Bachmann & Bergantz, 2004). For porous melt flow, it was also suggested that it might operate for millions of years: in the mid-crustal and lower-crustal migmatites in the Bohemian Massif, the duration of the porous melt flow was estimated to ~8–10 and ~15 Myr, respectively (U-Pb monazite geochronology; Hasalová et al., 2008c; Závada et al., 2021). In the Vosges Mts, a ~20 Myr-long melting period accompanied by porous melt flow in the middle crust was estimated (U-Pb monazite geochronology; Hasalová et al., 2015; Tabaud et al., 2014).
#### 4.3.5. Presence or Absence of Brittle Structures

Regarding brittle deformation, our models provide incomplete information as they simulate ductile (viscous) creep only. However, we can assess the potential for the opening of fractures using the compaction pressure that is one of the model variables. The calculated compaction pressure shows that fracturing is unlikely in the interior of the PMZ that is weakened by melt and deforms ductilely. Brittle deformation may occur mainly at the top margin of the PMZ/MZ where melt percolates into the relatively high-viscosity region above. This mechanism of brittle fracturing is more powerful if the rock permeability is high (see Figure 9a), which is the situation corresponding to meso-scale porous flow in contrast to the low-permeability micro-scale flow. The tendency to brittle deformation may be underestimated in our models because they neglect tectonic stresses and local stress build-up due to rock heterogeneity.

From the geological record, we see that melt-induced brittle or ductile microfracturing can occur in suprasolidus migmatites (Brown, 2004; Etheridge et al., 2021; Schulmann et al., 2008; Weinberg & Regenauer-Lieb, 2010; Závada et al., 2007). Within the rock, these mostly tensional microfractures occur predominantly in rheologically stronger mineral phases and for low melt fractions (e.g., Passchier & Trouw, 2005), and their distribution is controlled by the stress regime, melt pressure and rock anisotropy (e.g., Weinberg & Regenauer-Lieb, 2010). At a larger scale, the presence of dykes is mostly related to subsolidus conditions (e.g., Lister & Kerr, 1991). The lack of macro-scale brittle fractures is a typical feature of migmatitic terranes where micro-scale porous melt flow was described. In general, the occurrence of brittle structures at low melt fractions, at the margins of the partially molten zones and preferably in high-viscosity material is in agreement with the model results.

### 4.4. Model Limitations

The presented model is rather simple and in several respects it cannot capture the behavior of the continental crust. Below we discuss the most important limitations of our modeling approach: uncertainty in the determination of the Darcy coefficient and neglection of meso-scale melt segregation, simplified melting parametrization, absence of tectonic stresses, neglection of brittle deformation, inaccurate treatment of “magmatic” zones and issues related to the numerical formulation of the model.

In our models, we have identified the Darcy coefficient as the key parameter that determines the style of melt migration. This coefficient is the ratio of the rock permeability and melt viscosity, both of which are relatively weakly constrained and may vary significantly. We used the constant value of $10^4$ Pa s for the melt viscosity. This value is a reasonable estimate for hydrous felsic melts at pressures over ~8 GPa (e.g., Holtz et al., 2001). However, the melt viscosity varies with temperature and due to varying water content also with pressure. Holtz et al. (2001) calculated that the melt viscosity can vary by one order of magnitude over the pressure range of 0.5–0.8 GPa and temperature range of 700°C–950°C. They used two competing models for the viscosity, which resulted in similar average values and a clear decrease in the viscosity with the pressure. However, the way the viscosity varied with the temperature was quite different in the two models. Based on their results we expect that the melt viscosity will decrease with depth in the crust, but that this trend may be reversed in the coldest upper part of the PMZ due to the increased water content. While the change of one order of magnitude in the melt viscosity is relatively small from the experimentalists' point of view, it translates into a crucial change of style of melt migration and PMZ deformation in our models. A more accurate determination of the melt viscosity is therefore essential if we want to understand the melt migration processes in the hot continental crust.

The determination of the rock permeability is even more complicated. The porosity dependence of the permeability is not well constrained. Namely the term limiting the permeability for high porosities ($(1 - \phi)^2$ in our models) may be defined differently and in some studies it has been omitted entirely (e.g., Schmeling et al., 2019). We tested the effect of this term (Text S3 in Supporting Information S1) showing that it affects the formation of MZs but does not change the style of the model evolution. There is also a debate on the value of the permeability prefactor $k_0$. The choice of the permeability prefactor at the upper margin of the range of the reported values is justified by the observed geochemistry of mid-ocean basalts (e.g., Connolly et al., 2009). Still, the uncertainty of determination of this value is significant.

The assumption of a constant reference Darcy coefficient in each model requires a constant characteristic spacing of melt-filled pores, which contrasts with the spatially variable and dynamically changing melt distribution in many rocks. First, if the melt is distributed along grain boundaries, the change of the grain size, for example, due
to recrystallization affects the permeability (Equation 4). Second, melt segregation into a network of leucosomes or conduits would result in an increase of the characteristic spacing of the pores and therefore an increase of $K_{D0}$. This meso-scale segregation process cannot be resolved in our model and its physical description is beyond the scope of the present study. Based on the models we can therefore argue that melt can move through the crust by porous flow, but we cannot determine why in some cases it remains spread along grain boundaries, while in other cases it segregates into a meso-scale network. Third, in layered rocks, the anisotropy of the permeability may be important, because the melt flow in the direction perpendicular to the layering may be inhibited.

Keeping these limitations in mind, we argue that our models with a low-intermediate Darcy coefficient and low-moderate melt fractions describe well the behavior of rocks with a relatively homogeneous grain-scale network of pores, as that observed in metaigneous rocks with a low melt fraction. Melt segregation is hindered in these rocks and brittle fractures shall not focus melt either, because the effective viscosity in the interior of the partially molten region is low. However, brittle failure may occur at the contact of the partially molten region with the stronger unmolten crust. In these places, a network of fractures can develop that will have a higher permeability and therefore different melt-migration dynamics. Models with a high Darcy coefficient correspond to crust with a developed meso-scale network of pores; this approach was proposed by Schmeling et al. (2019). These models, however, represent only a crude approximation, because they neglect heterogeneity and anisotropy of the pore distribution that is typically larger in rocks with meso-scale melt zones than in rocks with pores along grain boundaries.

Another important limitation of the model is in the simplified melting parametrization. This linear melting law may lead to unrealistic compositions, which is visible mainly in models with high $K_{D0}$ and large $\Delta T_{SL}$ (e.g., models KD-10, DTSL500Ta550). In these models, enriched melt percolates into the pristine material, changing the bulk composition toward enrichment, which in turn promotes melting, further melt enrichment and melt propagation toward the surface. This feedback is partly artificial because, in our parametrization, the solidus and liquidus curves are parallel and do not intersect at the minimum (end-member) composition. For the same reason, also the interpretation of the model results in terms of rock composition (metasediments vs. metaigneous rocks) and its effect on melting and melt flow style should be taken with caution. A parametrization that would better describe the dependence of melting on composition, especially on the water content, shall be the topic of a future study. In addition, a melt composition that would reflect its water content could be used to parametrize the viscosity of the melt.

Plate-tectonic stresses are not included in our models. Forced deformation would inevitably affect the partially molten crust, and the interplay between the localization of deformation and melt transfer has to be explored if we want to compare the model results with real rock structures. In addition, deformation coupled with erosion may bring the hot lower crustal material to shallower levels, as commonly recorded in orogens, leading to different pressure-temperature paths recorded by the rocks. Tectonic forcing may also induce brittle failure in the cold, high-viscosity crust. The brittle faults and shear zones may then play the role of melt conduits through this portion of the crust. Despite significant advances in this field (Keller et al., 2013), accurate treatment of brittle deformation together with melt percolation remains a great challenge for future numerical models. In this study, we are limited by the simplified (purely viscous) rheology. This only allows us to identify regions with the potential for brittle failure resulting from high melt overpressure.

Our models are intended to illustrate the character of melt migration in the partially molten crust, but they are not targeted to study the formation and dynamics of MZs. For several reasons, these regions are not described realistically in the present model setup. First, the adopted mathematical description of two-phase flow is appropriate for a solid matrix with melt-filled pores, not for disaggregated solid carried by melt. Second, the lower cut-off value of shear viscosity (10¹⁶ Pa s) is several orders of magnitude higher than the real viscosity of magma (see e.g., Sparks et al., 2019). Third, much finer time and spatial resolutions are required to model the internal dynamics of MZs. Therefore, we do not interpret the internal structure of regions with the melt fraction above ~30%.

The numerical accuracy of the two-phase flow models in ASPECT was previously tested and reported on several model setups (Dannberg et al., 2019; Dannberg & Heister, 2016). We tested the newly implemented treatment of composition on simple cases such as melting of a compositionally anomalous region with/without melt flow. We tested the effect of resolution on the modeling results (Text S3 and Figure S1 in Supporting Information S1). In higher-resolution models, the character of the flow remains the same, but some features, such as the PWs and channels, are smaller (narrower) and focus more deformation and melt flow. A similar phenomenon was described by Dannberg et al. (2019) in the case of melt focusing toward an oceanic ridge.
With a lower resolution, numerical oscillations of porosity and other properties may occur. Based on these tests, we find the reference resolution sufficient because the style of model evolution in the reference and high-resolution models is the same. However, we should note that the resolution of the PWs and channels in our models is not optimal. In order to accurately model the PWs, one needs to resolve the compaction length (e.g., Dohmen & Schmeling, 2021). With the finest grid spacing of 390 m, the resolution with second-order elements is ~200 m, while the typical compaction length in the PMZ in the reference model is ~100–300 m. The shape and size of PWs in such a case cannot be represented accurately, as indeed confirmed by our tests. On the other hand, in the models with a high Darcy coefficient (model KD-10), the compaction length is ~500–1,000 m and the resolution is sufficient. Comparing the PWs in model KD-10 and in the reference model, we observe that they look very similar except for their larger sizes in model KD-10. This gives us confidence that the observed PWs are not mere numerical artifacts of poor spatial resolution. On the other hand, no PWs develop in the models with a low Darcy coefficient (model KD-12), but we can not exclude their potential formation at a sub-resolution scale.

Our tests further demonstrate the role of the minimum cut-off value of viscosity and of selected material parameters ($\alpha_c$, $\alpha_\phi$, $\Delta\rho_c$, $\log(1/\phi)$) instead of $1/\phi$ dependence of the bulk viscosity and porosity dependence of $k_\phi$ (Figures S2–S7). The minimum cut-off value of viscosity is $10^{16}$ Pa s in the main series of models. In our tests, if this parameter exceeds the typical viscosity in the PMZ (~$10^{17}$ Pa s), the sizes and shapes of the PWs are affected and the development of the MZ is largely exaggerated (Figure S2). The tests of the other material parameters show that their change may somewhat shift the boundaries between the different styles of evolution, but will not affect the general conclusions of our study. The most significant effect was observed in the case of the compositional dependence of density ($\Delta\rho_c$). With a reduced $\Delta\rho_c$, depleted (or less enriched) material can rise more easily due to thermal and melt-related buoyancy, and convective regime is preferred.

### 4.5. Comparison With Previous Modeling Studies of Melt Migration and Melting in Crust

Several numerical-modeling studies attempted to model melt-migration dynamics in the continental crust. So far the most advanced study of the interplay between the melt flow and solid deformation was performed by Keller et al. (2013). Their 2D models of melt transport include visco-elasto-plastic rock rheology and capture different styles of melt extraction: diapirism, visco-plastic decompaction channels and elasto-plastic dykes. However, their models lack melting/solidification and thermal evolution, both of which are essential for our study.

The interplay between the composition, temperature, melting and melt flow was studied, for example, by Jackson et al. (2005) and Solano et al. (2014) using 1D models. The compositional variations, which they reported for the column of material heated from below, are in broad agreement with our results. Namely, both models show that compositional equilibration between the migrating melt and the solid rock leads to a higher melt fraction in cooler shallower levels than in the hot but depleted lower crust. In contrast to their models, ours include lateral movements that complicate the simple 1D image and we observe that the upwelling hot partially molten material can form a diapir or induce melt-enhanced convection. Analysis of the solved equations shows that the 1D models of porous flow are valid only for a thin partially molten region because the potential for the development of convection depends on the third power of the vertical extent of the region (e.g., Schmeling et al., 2019). For example, in the reference model, the melt fraction of 10% and the viscosity of $10^{18}$ Pa s, the critical thickness for the development of convection is ~10 km. This thickness may be reached easily in thick and hot orogens.

Convection in partially molten crust was studied by Babeyko et al. (2002) who argued that it may be a viable mechanism to warm up the thickened continental crust, in contrast to pure melt migration without solid mass flow. This is confirmed by our study which shows a significantly warmer middle crust in models where solid material participates in the motion than in models where melt moves through a (almost) static crust. In contrast to Babeyko et al., we take into account the depletion/enrichment of the material and its effect on the melting temperature. As a result, the position of the PMZ moves from the lower crust toward the middle crust, because its vertical extent is given by the solidus of the depleted material in the lower crust and that of the enriched material in the middle crust (Leitch & Weinberg, 2002). Due to the temperature difference between the lower and middle crust, the total volume of melt progressively decreases. Therefore, convection is only a transient stage and there is no (quasi-) stationary convective state that could be reached.

The general setup of our models is similar to the models of Schmeling et al. (2019) but there are several important differences. In their models, the elevated thermal regime is boosted by an initial temperature anom-
ally in the middle-lower crust. This results in high melt fractions at the beginning of the simulation and quick formation of a MZ in the middle crust. In our case, the model evolution begins with a normal temperature profile and then we impose a long-lasting high-temperature anomaly at the bottom boundary. As a result, there is less melt and the MZ is smaller or absent. Schmeling et al. (2019) described the formation of two enrichment layers in the middle crust. These two layers correspond to the MZ and the top of the PMZ in our model.

In their recent study, Schmeling et al. (2023) assumed a more moderate thermal regime: they started their models from a conductive temperature profile and prescribed several heat pulses at the bottom boundary. With these conditions, they obtained two main modes of model evolution, which they called “batholith emplacement” and “convective recycling” modes. The batholith emplacement mode developed preferably in models where the heating took the form of several heat pulses, for a high Darcy coefficient and a high solid viscosity. In this mode, the relative velocity between the melt and the solid is higher and melt-rich caps (“batholiths”) form and persist in the middle crust on top of larger partially molten regions. In contrast, the convective recycling mode likely develops due to a single (but lower in amplitude) heat pulse, for a low Darcy coefficient and a low solid viscosity. It shows lower melt fractions (<20%) without melt-rich caps and with a significant component of solid velocity. Our results agree well with these results: in our models, PWs through a mostly static solid develop for a high Darcy coefficient and a high solid viscosity, just in opposite to the melt-enhanced convection. In addition, Schmeling et al. (2023) showed that the batholith emplacement mode is supported by a significant density difference between the reference and enriched compositions, in agreement with our results.

An important difference between our results and those of Schmeling et al. (2023, 2019) is in the absence of PWs in the latter ones. This might be the effect of insufficient resolution, but in all these models, the resolution is probably close to the limit that allows for the formation of PWs. Other possible reasons for the lack of PWs are in the higher viscosity cut-off value, or in a different porosity dependence of the solid viscosity and permeability. The exact definition of these material parameters strongly affects melt migration and formation of PWs. The impact of different porosity dependencies of these material parameters on the evolution of the partially molten zones would deserve a separate study.

## 5. Conclusions

Our numerical models suggest that porous melt flow can be a viable and long-lasting mechanism of melt transport in continental crust. Assuming realistic material properties, the melt can migrate in hot and thick continental crust through pores with a characteristic spacing of ~1 mm and larger. Despite its low rate, such porous flow can gradually form a large partially molten zone in the middle and lower crust and thus have a profound effect on their thermal state, deformation and composition.

Our results are consistent with the concept of large and long-lasting partially molten (mush) zones in continental crust. The models confirm the tight relationship between the partially molten zones and porous melt flow. They show that MZs may form on top of the partially molten zones, but they are not necessarily prominent or extensive.

The models support the observation that the activity of micro-scale porous flow over large distances is more likely in metaigneous rocks compared to metasediments. This is thanks to the slowly increasing melting curves typical for metaigneous rocks which, together with their lack of a strong anisotropy, impede melt segregation and the development of meso-macro-scale melt conduits.

The character of the simulated melt flow and the partially molten zones strongly depends on the model parameters. We examined the role of the rock permeability, melt and solid viscosities, melting parameters, and temperature conditions. Depending on these parameters, we obtained various styles of melt distribution, melt flow and solid deformation, and the related thermal and compositional evolution. We categorized these styles as melt-enhanced convection, growth of partially molten diapirs and melt percolation in PWs.

The present study provides a new perspective on melt transport in hot continental crust that can be compared with and tested against geological and geophysical data. We believe that it will stimulate new ideas and questions that will improve current conceptual models of crustal melting.
## Data Availability Statement

Open-source code ASPECT version 2.3.0 was used for all the finite element computations. Additional code and input scripts required for reproduction of the presented results are available in Maierová et al. (2023). Any other data were not used, nor created for this research.

## Acknowledgments

We thank Juliane Dannberg, Klára Kalousová, and Ondřej Šrámek for discussions on the modeling of two-phase flow. We are grateful to Harro Schmeling and an anonymous reviewer for their insightful and detailed reviews, and to Roberto Weinberg for careful reading and commenting on the manuscript. Their input greatly improved the quality of the paper. This research was supported by Grant 23-07821S of the Czech Science Foundation. All simulations were performed using the ASPECT software. Paraview and Inkscape were used to visualize the results.

## References

Aguilar, C., Štípská, P., Chopin, F., Schulmann, K., Pitra, P., Závada, P., et al. (2020). Syn-deformational melt percolation through a high-pressure orthogneiss and the exhumation of a subducted continental wedge (Orlica-Šniežnik Dome, NE Bohemian Massif). *International Journal of Earth Sciences*, 109(4), 1213–1246.

Annen, C., Blundy, J. D., Leuthold, J., & Sparks, R. S. J. (2015). Construction and evolution of igneous bodies: Towards an integrated perspective of crustal magmatism. *Lithos*, 230, 206–221.

Babeyko, A. Y., Sobolev, S. V., Trumbull, R., Oncken, O., & Lavier, L. (2002). Numerical models of crustal scale convection and partial melting beneath the Altiplano-Puna plateau. *Earth and Planetary Science Letters*, 199(3–4), 373–388.

Bachmann, O., & Bergantz, G. W. (2004). On the origin of crystal-poor rhyolites: Extracted from batholithic crystal mushes. *Journal of Petrology*, 45(8), 1565–1582.

Bangerth, W., Dannberg, J., Fraters, M., Gassmöller, R., Glerum, A., Heister, T., et al. (2020). ASPECT: Advanced solver for problems in Earth's ConvecTion. *User Manual*. https://doi.org/10.6084/m9.fishshare.4865333

Bangerth, W., Dannberg, J., Gassmoeller, R., & Heister, T. (2020). Aspect v2.2.0. Zenodo. https://doi.org/10.5281/zenodo.3924604

Beck, S. L., & Zandt, G. (2002). The nature of orogenic crust in the central Andes. *Journal of Geophysical Research: Solid Earth*, 107(B10), ESE7-1-ESE7-16. https://doi.org/10.1029/2000JB000124

Brown, M. (2004). The mechanism of melt extraction from lower continental crust of orogens. *Earth and Environmental Science Transactions of the Royal Society of Edinburgh*, 95(1–2), 35–48.

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. *Canadian Journal of Earth Sciences*, 47(5), 655–694.

Cao, W., Kaus, B. J., & 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. https://doi.org/10.1002/2015TC004076

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

Cavalcante, C., Hollanda, M. H., Vauchez, A., & Kawata, M. (2018). How long can the middle crust remain partially molten during orogeny? *Geology*, 46(10), 839–842.

Christensen, N. I., & Mooney, W. D. (1995). Seismic velocity structure and composition of the continental crust: A global view. *Journal of Geophysical Research: Solid Earth*, 100(B6), 9761–9788.

Clark, M. K., & Royden, L. H. (2000). Topographic ooze: Building the eastern margin of Tibet by lower crustal flow. *Geology*, 28(8), 703–706.

Clemens, J. (2006). *Melting of the continental crust: Fluid regimes, melting reactions, and source-rock fertility*. Cambridge University Press.

Clemens, J., & Vielzeuf, D. (1987). Constraints on melting and magma production in the crust. *Earth and Planetary Science Letters*, 86(2–4), 287–306.

Connolly, J., & Podladchikov, Y. Y. (2007). Decompression weakening and channeling instability in ductile porous media: Implications for asthenospheric melt segregation. *Journal of Geophysical Research: Solid Earth*, 112(B10). https://doi.org/10.1029/2005JB004213

Connolly, J., Schmidt, M. W., Solferino, G., & Bagdassarov, N. (2009). Permeability of asthenospheric mantle and melt extraction rates at mid-ocean ridges. *Nature*, 462(7270), 209–212.

Daczko, N. R., Piazolo, S., Meek, U., Stuart, C. A., & Elliott, V. (2016). Hornblendite delineates zones of mass transfer through the lower crust. *Scientific Reports*, 6(1), 1–6.

Dannberg, J., Gassmöller, R., Grove, R., & Heister, T. (2019). A new formulation for coupled magma/mantle dynamics. *Geophysical Journal International*, 219(1), 94–107.

Dannberg, J., & Heister, T. (2016). Compressible magma/mantle dynamics: 3D, adaptive simulations in ASPECT. *Geophysical Journal International*, 207(3), 1343–1366. https://doi.org/10.1093/gji/ggw329

Delph, J. R., Ward, K. M., Zandt, G., Ducea, M. N., & Beck, S. L. (2017). Imaging a magma plumbing system from MASH zone to magma reservoir. *Earth and Planetary Science Letters*, 457, 313–324.

del Potro, R., Díez, M., Blundy, J., Camacho, A. G., & Gottsmann, J. (2013). Diapiric ascent of silicic magma beneath the Bolivian Altiplano. *Geophysical Research Letters*, 40(10), 2044–2048. https://doi.org/10.1002/grl.50493

Ding, H., Kohn, M. J., & Zhang, Z. (2021). Long-lived (ca. 22–24 Myr) partial melts in the eastern Himalaya: Petrochronologic constraints and tectonic implications. *Earth and Planetary Science Letters*, 558, 1167–1176.

Dohmen, J., & Schmeling, H. (2021). Magma ascent mechanisms in the transition regime from solitary porosity waves to diapirism. *Solid Earth*, 12(7), 1549–1561.

Dong, H., Wei, W., Jin, S., Ye, G., Jones, A. G., Zhang, L., et al. (2020). Shaping the surface deformation of central and south Tibetan Plateau: Insights from magnetotelluric array data. *Journal of Geophysical Research: Solid Earth*, 125(9), e2019JB019206. https://doi.org/10.1029/2019JB019206

Etheridge, M. A., Daczko, N. R., Chapman, T., & Stuart, C. A. (2021). Mechanisms of melt extraction during lower crustal partial melting. *Journal of Metamorphic Geology*, 39(1), 57–75.

Gaillard, F., Scaillet, B., & Pichavant, M. (2004). Evidence for present-day leucogranite pluton growth in Tibet. *Geology*, 32(9), 801–804.

Giordano, D., & Dingwell, D. B. (2003). Non-Arrhenian multicomponent melt viscosity: A model. *Earth and Planetary Science Letters*, 208(3–4), 337–349.

Hacker, B., Ritzwoller, M., & Xie, J. (2014). Partially melted, mica-bearing crust in Central Tibet. *Tectonics*, 33(7), 1408–1424. https://doi.org/10.1002/2014TC003545

Hall, D., & Kisters, A. (2012). The stabilization of self-organised leucogranite networks—Implications for melt segregation and far-field melt transfer in the continental crust. *Earth and Planetary Science Letters*, 355, 1–12.

Hasalová, P., Janoušek, V., Schulmann, K., Štípská, P., & Erban, V. (2008). From orthogneiss to migmatite: Geochemical assessment of the melt infiltration model in the Gföhl Unit (Moldanubian Zone, Bohemian Massif). *Lithos*, 102(3–4), 508–537.
Hasalová, P., Schulmann, K., Lexa, O., Štípská, P., Hrouda, F., Ulrich, S., et al. (2008a). Origin of migmatites by deformation-enhanced melt infiltration of orthogneiss: A new model based on quantitative microstructural analysis. *Journal of Metamorphic Geology*, 26(1), 29–53.

Hasalová, P., Schulmann, K., Tabaud, A. S., & Oliot, E. (2015). Microstructural evidences for mineralogical inheritance in partially molten rocks: Example from the Vosges Mts. *Bulletin de la Societe Geologique de France*, 186(2–3), 131–143.

Hasalová, P., Štípská, P., Powell, R., Schulmann, K., Janoušek, V., & Lexa, O. (2008b). Transforming mylonitic metagranite by open-system interactions during melt flow. *Journal of Metamorphic Geology*, 26(1), 55–80.

Hasalová, P., Weinberg, R., & Macrae, C. (2011). Microstructural evidence for magma confluence and reusage of magma pathways: Implications for magma hybridization, Karakoram Shear Zone in NW India. *Journal of Metamorphic Geology*, 29(8), 875–900.

Hawkesworth, C., & Kemp, A. (2006). The differentiation and rates of generation of the continental crust. *Chemical Geology*, 226(3–4), 134–143.

Heister, T., Dannberg, J., Gassmöller, R., & Bangerth, W. (2017). High accuracy mantle convection simulation through modern numerical methods. II: Realistic models and problems. *Geophysical Journal International*, 210(2), 833–851. https://doi.org/10.1093/gji/ggx195

Holtz, F., Johannes, W., Tamic, N., & Behrens, H. (2001). Maximum and minimum water contents of granitic melts generated in the crust: A reevaluation and implications. *Lithos*, 56(1), 1–14.

Hyndman, R. (2019). Mountain building orogeny in precollision hot backarcs: North American Cordillera, India-Tibet, and Grenville Province. *Journal of Geophysical Research: Solid Earth*, 124(2), 2057–2079. https://doi.org/10.1029/2018JB016697

Jackson, M., Gallagher, K., Petford, N., & Cheadle, M. (2005). Towards a coupled physical and chemical model for tonalite–trondhjemite–granodiorite magma formation. *Lithos*, 79(1–2), 43–60.

Jamieson, R. A., Unsworth, M. J., Harris, N. B., Rosenberg, C. L., & Schulmann, K. (2011). Crustal melting and the flow of mountains. *Elements*, 7(4), 253–260.

Kelemen, P. B., Dick, H. J., & Quick, J. E. (1992). Formation of harzburgite by pervasive melt/rock reaction in the upper mantle. *Nature*, 358(6388), 635–641.

Keller, T., May, D. A., & Kaus, B. J. (2013). Numerical modelling of magma dynamics coupled to tectonic deformation of lithosphere and crust. *Geophysical Journal International*, 195(3), 1406–1442.

Korhonen, F., Saito, S., Brown, M., & Siddoway, C. (2010). Modeling multiple melt loss events in the evolution of an active continental margin. *Lithos*, 116(3–4), 230–248.

Kronbichler, M., Heister, T., & Bangerth, W. (2012). High accuracy mantle convection simulation through modern numerical methods. *Geophysical Journal International*, 191, 12–29. https://doi.org/10.1111/j.1365-246X.2012.05609.x

le Breton, N., & Thompson, A. B. (1988). Fluid-absent (dehydration) melting of biotite in metapelites in the early stages of crustal anatexis. *Contributions to Mineralogy and Petrology*, 99(2), 226–237.

Leitch, A., & Weinberg, R. F. (2002). Modelling granite migration by mesoscale pervasive flow. *Earth and Planetary Science Letters*, 200(1–2), 131–146.

Lister, J. R., & Kerr, R. C. (1991). Fluid-mechanical models of crack propagation and their application to magma transport in dykes. *Journal of Geophysical Research: Solid Earth*, 96(B6), 10049–10077.

Lu, L. X., & Jiang, D. (2019). Quartz flow law revisited: The significance of pressure dependence of the activation enthalpy. *Journal of Geophysical Research: Solid Earth*, 124(1), 241–256. https://doi.org/10.1029/2018JB016226

Maierová, P., Hasalová, P., Schulmann, K., Štípská, P., & Souček, O. (2023). Porous melt flow in continental crust—A numerical modeling study [Software]. Zenodo. https://doi.org/10.5281/zenodo.8181433

McKenzie, D. (1984). The generation and compaction of partially molten rock. *Journal of Petrology*, 25(3), 713–765.

Meek, U., Piazolo, S., & Daczko, N. R. (2019). The field and microstructural signatures of deformation-assisted melt transfer: Insights from magmatic arc lower crust, New Zealand. *Journal of Metamorphic Geology*, 37(6), 795–821.

Passchier, C. W., & Trouw, R. A. (2005). *Microtectonics*. Springer Science & Business Media.

Patiño Douce, A. E., & Harris, N. (1998). Experimental constraints on Himalayan anatexis. *Journal of Petrology*, 39(4), 689–710.

Petford, N., Cruden, A., McCaffrey, K., & Vigneresse, J.-L. (2000). Granite magma formation, transport and emplacement in the Earth’s crust. *Nature*, 408(6813), 669–673.

Pritchard, M. E., & Gregg, P. M. (2016). Geophysical evidence for silicic crustal melt in the continents: Where, what kind, and how much? *Elements*, 12(2), 121–127.

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

Reichardt, H., & Weinberg, R. F. (2012). Hornblende chemistry in meta-and diatexites and its retention in the source of leucogranites: An example from the Karakoram Shear Zone, NW India. *Journal of Petrology*, 53(6), 1287–1318.

Rosenberg, C., & Handy, M. (2005). Experimental deformation of partially melted granite revisited: Implications for the continental crust. *Journal of Metamorphic Geology*, 23(1), 19–28.

Rudge, J. F. (2018). The viscosities of partially molten materials undergoing diffusion creep. *Journal of Geophysical Research: Solid Earth*, 123(12), 10–534. https://doi.org/10.1029/2018JB016530

Rudnick, R., & Gao, S. (2003). The composition of the continental crust. In H. D. Holland & K. K. Turekian (Eds.), *The crust: Treatise on geochemistry* (Vol. 3, pp. 1–64). Elsevier-Pergamon.

Rutter, E. H., Mecklenburgh, J., Brown, M., & Rushmer, T. (2006). *The extraction of melt from crustal protoliths and the flow behaviour of partially molten crustal rocks: An experimental perspective*. Citeseer.

Sandiford, M., & McLaren, S. (2002). Tectonic feedback and the ordering of heat producing elements within the continental lithosphere. *Earth and Planetary Science Letters*, 204(1–2), 133–150.

Sawyer, E. (2001). Melt segregation in the continental crust: Distribution and movement of melt in anatectic rocks. *Journal of Metamorphic Geology*, 19(3), 291–309.

Sawyer, E., Cesare, B., & Brown, M. (2011). When the continental crust melts. *Elements*, 7(4), 229–234.

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. *Geophysical Journal International*, 190(3), 1571–1578.

Schmeling, H., Marquart, G., Weinberg, R., & Kumaravel, P. (2023). Dynamic two-phase flow modelling of melt segregation in continental crust: Batholith emplacement vs. crustal convection, with implications for magmatism in thickened plateaus. *Geochemistry, Geophysics, Geosystems*, 24(5). https://doi.org/10.1029/2023GC010860

Schmeling, H., Marquart, G., Weinberg, R., & Wallner, H. (2019). Modelling melting and melt segregation by two-phase flow: New insights into the dynamics of magmatic systems in the continental crust. *Geophysical Journal International*, 217(1), 422–450.

Schulmann, K., Konopásek, J., Janoušek, V., Lexa, O., Lardeaux, J.-M., Edel, J.-B., et al. (2009). An Andean type Palaeozoic convergence in the Bohemian massif. *Comptes Rendus Geoscience*, 341(2–3), 266–286.
Schulmann, K., Martelat, J.-E., Ulrich, S., Lexa, O., Štípská, P., & Becker, J. K. (2008). Evolution of microstructure and melt topology in partially molten granitic mylonite: Implications for rheology of felsic middle crust. *Journal of Geophysical Research: Solid Earth*, 113(B10). https://doi.org/10.1029/2007JB005508

Schwindinger, M., & Weinberg, R. F. (2017). A felsic MASH zone of crustal magmas—Feedback between granite magma intrusion and in situ crustal anatexis. *Lithos*, 284, 109–121.

Simpson, G., Spiegelman, M., & Weinstein, M. I. (2010). A multiscale model of partial melts: 2. Numerical results. *Journal of Geophysical Research: Solid Earth*, 115(B4). https://doi.org/10.1029/2009JB006376

Solano, J., Jackson, M., Sparks, R., & Blundy, J. (2014). Evolution of major and trace element composition during melt migration through crystalline mush: Implications for chemical differentiation in the crust. *American Journal of Science*, 314(5), 895–939.

Sparks, R., Annen, C., Blundy, J., Cashman, K., Rust, A., & Jackson, M. (2019). Formation and dynamics of magma reservoirs. *Philosophical Transactions of the Royal society A*, 377(2139), 20180019.

Spiegelman, M. (1993). Physics of melt extraction: Theory, implications and applications. *Philosophical Transactions of the Royal Society of London, Series A: Physical and Engineering Sciences*, 342(1663), 23–41.

Štípská, P., Hasalová, P., Powell, R., Závada, P., Schulmann, K., Racek, M., et al. (2019). The effect of melt infiltration on metagranitic rocks: The Snieznik dome, Bohemian Massif. *Journal of Petrology*, 60(3), 591–618.

Stuart, C. A., Piazolo, S., & Daczko, N. R. (2016). Mass transfer in the lower crust: Evidence for incipient melt assisted flow along grain boundaries in the deep arc granulites of Fiordland, New Zealand. *Geochemistry, Geophysics, Geosystems*, 17(9), 3733–3753. https://doi.org/10.1002/2015GC006236

Stuart, C. A., Piazolo, S., & Daczko, N. R. (2018). The recognition of former melt flux through high-strain zones. *Journal of Metamorphic Geology*, 36(8), 1049–1069.

Tabaud, A.-S., Whitechurch, H., Rossi, P., Schulmann, K., Guerrot, C., & Cocherie, A. (2014). Devonian–Permian magmatic pulses in the northern Vosges Mountains (NE France): Result of continuous subduction of the Rhenohercynian Ocean and Avalonian passive margin. *Geological Society, London, Special Publications*, 405(1), 197–223.

Tanner, D. C. (1999). The scale-invariant nature of migmatite from the Oberpfalz, NE Bavaria and its significance for melt transport. *Tectonophysics*, 302(3–4), 297–305.

Unsworth, M. (2010). Magnetotelluric studies of active continent–continent collisions. *Surveys in Geophysics*, 31(2), 137–161.

Vanderhaeghe, O. (1999). Pervasive melt migration from migmatites to leucogranite in the Shuswap metamorphic core complex, Canada: Control of regional deformation. *Tectonophysics*, 312(1), 35–55.

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

Vernon, R., & Paterson, S. (2001). Axial-surface leucosomes in anatectic migmatites. *Tectonophysics*, 335(1–2), 183–192.

Ward, K. M., Zandt, G., Beck, S. L., Christensen, D. H., & McFarlin, H. (2014). Seismic imaging of the magmatic underpinnings beneath the Altiplano-Puna volcanic complex from the joint inversion of surface wave dispersion and receiver functions. *Earth and Planetary Science Letters*, 404, 43–53.

Wark, D. A., & Watson, E. B. (1998). Grain-scale permeabilities of texturally equilibrated, monomineralic rocks. *Earth and Planetary Science Letters*, 164(3–4), 591–605.

Wei, W., Unsworth, M., Jones, A., Booker, J., Tan, H., Nelson, D., et al. (2001). Detection of widespread fluids in the Tibetan crust by magnetotelluric studies. *Science*, 292(5517), 716–719.

Weinberg, R. F., & Hasalová, P. (2015). Water-fluxed melting of the continental crust: A review. *Lithos*, 212, 158–188.

Weinberg, R. F., & Regenauer-Lieb, K. (2010). Ductile fractures and magma migration from source. *Geology*, 38(4), 363–366.

Weinberg, R. F., Vernon, R. H., & Schmeling, H. (2021). Processes in mushes and their role in the differentiation of granitic rocks. *Earth-Science Reviews*, 220, 103665.

Wen, D.-P., Wang, Y.-F., Zhang, J.-F., Li, P.-X., & Jin, Z.-M. (2021). Rheology of felsic granulite at high temperature and high pressure. *Journal of Geophysical Research: Solid Earth*, 126(2), e2020JB020966. https://doi.org/10.1029/2020JB020966

White, R., Powell, R., & Holland, T. (2001). Calculation of partial melting equilibria in the system Na₂O-CaO-K₂O-FeO-MgO-Al₂O₃-SiO₂-H₂O (NCKFMASH). *Journal of Metamorphic Geology*, 19(2), 139–153.

Whittington, A., Richet, P., Behrens, H., Holtz, F., & Scaillet, B. (2004). Experimental temperature–X (H₂O)–viscosity relationship for leucogranites and comparison with synthetic silicic liquids. *Earth and Environmental Science Transactions of the Royal Society of Edinburgh*, 95(1–2), 59–71.

Wolfram, L. C., Weinberg, R. F., Nebel, O., Hamza, K., Hasalová, P., Miková, J., & Becchio, R. (2019). A 60-Myr record of continental back-arc differentiation through cyclic melting. *Nature Geoscience*, 12(3), 215–219.

Yakymchuk, C., Brown, M., Ivanic, T., & Korhonen, F. (2013). Leucosome distribution in migmatitic paragneisses and orthogneisses: A record of self-organized melt migration and entrapment in a heterogeneous partially-molten crust. *Tectonophysics*, 603, 136–154.

Závada, P., Schulmann, K., Konopásek, J., Ulrich, S., & Lexa, O. (2007). Extreme ductility of feldspar aggregates—Melt-enhanced grain boundary sliding and creep failure: Rheological implications for felsic lower crust. *Journal of Geophysical Research: Solid Earth*, 112(B10). https://doi.org/10.1029/2006JB004820

Závada, P., Schulmann, K., Racek, M., Hasalová, P., Jeřábek, P., Weinberg, R. F., et al. (2018). Role of strain localization and melt flow on exhumation of deeply subducted continental crust. *Lithosphere*, 10(2), 217–238.

Závada, P., Štípská, P., Hasalová, P., Racek, M., Jeřábek, P., Schulmann, K., et al. (2021). Monazite geochronology in melt-percolated UHP meta-granitoids: An example from the Erzgebirge continental subduction wedge, Bohemian Massif. *Chemical Geology*, 559, 119919.

Zhang, Z., Deng, Y., Teng, J., Wang, C., Gao, R., Chen, Y., & Fan, W. (2011). An overview of the crustal structure of the Tibetan plateau after 35 years of deep seismic soundings. *Journal of Asian Earth Sciences*, 40(4), 977–989.