A continuum model of multi-phase reactive transports in igneous systems

Multi-phase reactive transports are ubiquitous in igneous systems. Complex mechanical and thermodynamic interactions between solid, liquid and gas phases in magmatic aggregates lead to strongly non-linear processes including flow localization, reaction-transport feedbacks, and tipping-point behaviors. A challenging aspect of igneous processes is that they range from solid- dominated porous to liquid-dominated suspension flows and therefore entail a wide spectrum of rheological conditions, flow speeds, and length scales. Previous models have mostly been limited to the limits of two-phase flow in competent rock at low melt fraction and in largely liquid magma bearing a small crystal load. The goal of this paper is to develop a framework that is able to capture all stages of an igneous system from source to surface including not only a rock and a melt but also an exsolved volatile phase. Adding a third phase is critical for the study of magma degassing, which is intimately linked with magma transport and storage in the shallow crust. Here, we derive an n-phase reactive transport model building on the general concepts of Mixture Theory and the specific principles of Rational Thermodynamics and Non-equilibrium Thermodynamics. The model operates at the continuum scale and requires constitutive relations for fluxes within and transfers between phases. While our description of fluxes is similar to previous work, we introduce a new formulation for phase transfers of mass, entropy, momentum, and volume. We further propose phenomenological closures for the material response coefficients that determine how fluxes and transfers respond to forcing gradients and phase differences in thermodynamic properties. Finally, we demonstrate that the known limits of two-phase porous and suspension flow emerge as special cases of our general model and discuss ramifications for igneous systems.


Introduction
The complex interplay between mechanics and thermodynamics in reactive transport processes is a common theme in many natural multi-phase systems. The intricate interactions of multiple material phases-solids, liquids, and gases-exhibit a wide range of processes pertinent to our fundamental understanding of the Earth, as well as to many applied science and engineering problems. Examples include methane seepage through thawing permafrost or along geological faults (e.g., Christensen et al., 2004;Etiope, 2009), the evolution of hydrocarbon reservoirs (e.g. Roure et al., 2005), and hydrothermal ore formation (e.g., Sillitoe, 2003Sillitoe, , 2010, as well as the related engineering problems of carbon sequestration (e.g. Gaus et al., 2005), hydraulic fracturing (e.g. Jha and Juanes, 2014), and ore concentration (e.g. Cariaga et al., 2005). Reactive multi-phase systems are highly nonlinear with reactions driving transport and vice versa. The result are complex feedbacks that may lead to localization of natural phenomena in space and time and chemical differentiation or selective element concentration.
The primary objective of this study is to develop a model that captures the key characteristics of reactive multi-phase flows in volcanic systems and their deep magmatic roots, collectively known as igneous systems. Improving our capacity to model igneous processes will advance our understanding of volcanic activity and associated hazards, of planetary differentiation that created a habitable Earth and the deep volatile cycles maintaining it (Lenardic et al., 2016, e.g.), and of the magmatichydrothermal generation of economic deposits of iron, copper, gold and other metals (Sillitoe, 2003(Sillitoe, , 2010. Despite focusing on igneous systems, we cast the model in a general form to allow application to other natural and engineering contexts more amenable to direct observations than the igneous systems buried deep in the subsurface. Developing a macroscopic model of coupled thermodynamic and mechanical behavior at scales much larger than microscopic phase constituents-solid grains, liquid films, tubes ore droplets, gas bubbles-builds on a long line of scientific inquiry (e.g., De Groot and Mazur, 1984;Truesdell, 1984). The most common approach for formulating multi-phase models is known as Mixture Theory. The theory assumes that macroscopic effects of microscopic phase interactions may be modeled as an interacting set of interpenetrating continuum fields that represent the mixture of phases at the system scale. Implicit to the approach is the concept of scale separation, namely that it is possible-indeed, appropriate-to spatially average over processes at the local scale to derive a simplified model of system-scale behavior without capturing scales in between (e.g., Anderson and Jackson, 1967;Slattery, 1967).
In mixture models, the distributed microscopic properties of a representative volume of multiphase aggregate are represented by a point on a set of continuum fields (see Fig. 1). On the local scale, phase constituents including solid grains, liquid films, and gas bubbles occupy a finite volume and interact across well-defined interfaces. On the continuum scale, they are represented by the volume fraction they occupy in a control volume and their behavior is described by averaged phase interaction terms. Needless to say, different phase distributions may result in identical phase fractions, and system-scale interaction terms may not fully capture local processes, highlighting that some information will be lost in the transition to the continuum scale.
Formally, the transition of scales is achieved by spatial averaging. Multiple techniques have been proposed to accomplish the task (e.g., Boruvka et al., 1985;Quintard and Whitaker, 1988;Gray et al., 1993;Whitaker, 2013). The aim of these methods is to replace the complexity of local-scale topologies and processes with appropriately defined spatially averaged terms. For example, spatial averaging would be used to represent the fluctuating motion of phase constituents within a control volume by a single averaged phase velocity. Building a mixture model requires the formulation of conservation statements and thermody-  Figure 1: Schematic of a control volume (cube) with phase constituents in an igneous aggregate at the local scale and corresponding points (circles) on set of continuum phase fields (squares) at the system scale. Image is a microtomographic scan of a partially molten rock modified from Zhu et al. (2011); vapour bubbles (black) added to illustrate a three-phase system. Dashed and solid arrows denote fluxes within and transfers between phases. The latter represent system-scale effects of local-scale phase interactions, capture multi-phase behavior at continuum scale. namic principles, along with constitutive relations for system-scale processes, and closures prescribing how processes respond to applied forces as a function of material properties. Constitutive relations are required specifically for averaged transports of thermodynamic properties within phases between control volumes (fluxes) and between phases within a control volume (transfers) (see Fig. 1) as a function of system variables (e.g., phase velocities, pressures, temperatures, etc.). The material closures encapsulate how system-scale fluxes and transfers depend on local-scale material properties (e.g., density, viscosity, compressibility, etc.) and phase topologies (e.g., granular size, phase connectivity, wetting angles, etc.). Choosing these closures generally constitutes the most consequential step of model building in a continuum framework.
The main difference between mixture model approaches is the scale at which conservation statements, thermodynamic principles, constitutive relations, and material closures are formulated. In the Method of Volume Averaging (e.g., Whitaker, 1988, 1996;Whitaker, 2013), conservation equations are both formulated and closed at the local scale. The closed equations are then averaged to the system scale through a hierarchy of averages (Quintard and Whitaker, 1988;Whitaker, 2013). While powerful when applicable, this approach assumes a temporally fixed unit cell and is hence difficult to generalise to dynamic problems including suspension flows or porous flow in a deforming matrix.
An alternative approach is Thermodynamically Constrained Averaging Theory (e.g., Gray et al., 2013;Gray and Miller, 2014;Dye et al., 2015), which writes the governing equations at the con-tinuum scale but formulates constitutive relations for mechanical and thermodynamic processes at the microscale. The advantage is that local-scale energy balances are derived according to first principles. This strategy is particularly powerful for complex wetting problems that depend sensitively on interface energetics (e.g., Miller, 2010, 2011a). The resulting thermodynamic constraints are then used to derive local-scale equilibrium conditions (Boruvka et al., 1985) that inform the averaging procedure to the system scale (Gray et al., 1993;Miller and Gray, 2005). As a consequence, system-scale variables are formally linked to local-scale variables.
Thermodynamically Constrained Averaging perhaps comes closest to providing a rigorous and general multi-scale model framework for reactive multi-phase transport, at least in the porous flow limit where it has mostly been applied (e.g., Gray and Miller, 2014;Dye et al., 2015). The resulting equations, however, become somewhat unwieldy. In addition to the standard conservation laws, they entail evolution equations for variables related to local-scale phase topology including common curves and triple points (e.g., Gray and Miller, 2010), as well as interfacial forces including surface tensions and capillary pressures (e.g., Gray and Miller, 2011b). In igneous systems, detailed observations of microscale dynamics are oftentimes difficult to obtain and, if available, are complex and challenging to interpret. Hence, for a leading-order physical understanding, the additional mathematical complexity required to formally couple the scales may not always be justified.
We opt to follow a third approach, inspired by the practical success of two particular mixture models: the two-phase porous flow model of McKenzie (1984) and the two-fluid model of Drew (1971). Both follow the Rational Thermodynamics approach (e.g., Truesdell, 1984;Drew and Passman, 2006), which averages local-scale conservation equations for phases and interfaces to the system scale before they are closed (e.g., Hassanizadeh and Gray, 1979;Lahey and Drew, 1988). Thermodynamic principles are then formulated at the continuum scale and exploited to constrain admissible constitutive relations (e.g., Müller, 1967Müller, , 1968Liu, 1972;Hassanizadeh and Gray, 1979). Hence, this class of models operates almost exclusively at the continuum scale. The key drawback is that it forgoes a first-principles link to local-scale variables and processes (e.g., Gray, 1980, 1990;Svendsen and Hutter, 1995;Bennethum et al., 2000). While the approach yields expression resembling familiar single-phase equations, the disconnect between well-defined thermodynamic variables at the local scale and their more loosely defined system-scale counterparts entails terms that are not straightforward to constrain or interpret. The system-scale constitutive relations and material closures employed in these models, while loosely constrained by thermodynamic principles, are inevitably non-unique and phenomenological in nature.
In summary, our choice of following the continuum-scale approach is motivated by the pragmatic aim of providing a mathematically tractable and intuitively accessible framework for building phenomenological yet physically and thermodynamically consistent models. In the following, we start our derivation from a generic conservation equation at the continuum scale, which we use to formulate conservatoin principles for mass, momentum and energy. We then use thermodynamic principles to assemble an expanded entropy inequality, which we exploit to choose admissible constitutive relations. This procedure was pioneered in Non-equilibrium Thermodynamics (De Groot and Mazur, 1984;Jou et al., 2001) and formalised by Rational Thermodynamics (e.g., Truesdell, 1984;Drew and Passman, 2006). After proposing phenomenological closures for the material response coefficients arising in the chosen constitutive relations, we finally demonstrate how our model is applied to reactive transports from porous to suspension flows in igneous aggregates of two, three, or more phases.

Multi-phase reactive transports in igneous systems
Igneous processes occupy a broad range of temporal and spatial scales (Crisp, 1984;Caricchi, 2014;Papale, 2018). As non-linear systems, they can exhibit dramatic and sudden shifts in behavior. One consequential expressions of such tipping-point behavior occurs when silicic crustal magma reservoirs that have been stagnant for centuries to millennia are mobilised to renewed explosive activity within few years, weeks or even days (e.g. Burgisser and Bergantz, 2011;Cooper and Kent, 2014).
A uniquely challenging aspect of igneous systems is their wide range of mechanical regimes, from porous flow in the predominantly solid upper mantle to suspension flow in melt-rich magma bodies feeding into volcanic vents. Adding to the challenge is the complex thermodynamics of petrological phase equilibria spanning a wide compositional space including silicates, metal oxides, and volatile species (e.g., Cashman et al., 2017). To motivate our model framework, we briefly discuss arc magmatism as an example igneous system (see Fig. 2). We distinguish four stages that are pertinent to other igneous systems as well: (1) melt generation in the asthenosphere, (2) melt focusing into the lithosphere, (3) magma processing in the crust, and (4) magma degassing from the shallow crust to the surface. Hydrous fluids sourced from breakdown reactions of hydrated minerals percolate from the subducting plate (Cagnioncle et al., 2007;Wilson et al., 2014) and react with the asthenospheric mantle wedge to produce primitive arc melts (Gaetani and Grove, 1998;Grove et al., 2006). These melts segregate from their mantle residue by reactive porous flow at low melt fractions (Fig. 2, stage 1). McKenzie's model (McKenzie, 1984) and similar formulations (Fowler, 1985;Ribe, 1985a;Scott and Stevenson, 1986) based on earlier works including Drew (1971), Sleep (1974), and Turcotte and Ahern (1978) have been successfully employed to describe the porous flow of silicate melts through a compacting rock matrix. Bercovici et al. (2001) and Bercovici and Ricard (2003) introduced a phase-symmetrical generalization including surface tensions and a more rigorous treatment of the energetics of the problem. When coupled to models of thermo-chemical evolution and meltrock reactions (Šrámek et al., 2007;Rudge et al., 2011), non-linear reaction-transport feedbacks emerge such as the reactive infiltration instability (Chadam et al., 1986) leading to channelised melt transport (Aharonov et al., 1995;Spiegelman et al., 2001;Keller and Katz, 2016) that can have important ramifications for the geochemistry of melts and their solid residue (Spiegelman and Elliott, 1993;Spiegelman and Kelemen, 2003;Hewitt, 2010).
As mantle melts percolate into the cooler mantle lithosphere and crust, thermal exposure drives crystallization reactions that deposit latent heat, volatiles and fertile minerals into the mantle lithosphere (e.g., Keller et al., 2017) (Fig. 2, stage 2). Locations of high melt flux may develop hot funnels into which melts from the surrounding area are focused (England and Katz, 2010;Rees Jones et al., 2018). Throughout the crust, magma transport may stall along rheological and density contrasts (e.g., Ritter et al., 2013) and melt content may come to exceed the disaggregation threshold of the granular matrix (e.g. Arzi, 1978;Van der Molen and Paterson, 1979;Renner et al., 2000). Depending on the mechanical and thermal structure of the crust, its thickness (Hildreth and Moorbath, 1988), and regional tectonic regime (Cembrano and Lara, 2009), transport towards the shallow crust is thought to be accommodated by diapirs rising through ductile crust (Bateman, 1984;Cruden, 1990) or by fractures propagating through brittle rock (Clemens and Mawer, 1992;Rubin, 1993;Rozhko et al., 2007;Havlin et al., 2013;Keller et al., 2013).
Repeated injection of melt into the lower crust combined with crystallization, fractionation, and amalgamation of magma batches has been invoked (Hildreth and Moorbath, 1988;Annen et al., 2006;Cashman and Blundy, 2013) to explain differentiation from primitive to evolved magmas. The latter make up the bulk of eruptive products at major arc volcanic centers and form large plutonic bodies that generate continental crust. Fractional crystallization, i.e., the settling out of crystallisation products from melt-rich crustal magma chambers (Bowen, 1915) has long been the main paradigm of magmatic differentiation and has been investigated in a range of igneous suspension models (Huppert and Sparks, 1981;Brandeis and Jaupart, 1986;Martin and Nokes, 1988;Rudman, 1992;Bergantz and Ni, 1999;Dufek and Bachmann, 2010;Molina et al., 2012). Today, consensus is building towards a new paradigm of magma processing in trans-crustal crystalrich mush bodies with interspersed and likely transient crystal-poor magma lenses (Caricchi and Blundy, 2015;Cashman et al., 2017) (Fig. 2, stage 3). To test this hypothesis, models are required to bridge the limits of porous and suspension flows including the mush regime in between.
Finally, in the shallow crust, volatile exsolution from the melt phase marks the final stage of the igneous system: the degassing of sub-volcanic magma bodies that drives volcanic eruptions and may generate magmatic-hydrothermal ore deposits (Fig. 2, stage 4). Volatile exsolution upon magma decompression and cooling (Papale, 1999) produces droplets of supercritical fluids, brines, and/or bubbles of gaseous vapour (Bachmann and Bergantz, 2006;Driesner and Heinrich, 2007;Driesner, 2007;Ruprecht et al., 2008). However, volatile depletion of the melt may also increase crystallization (e.g., Métrich et al., 2001). Thus, magma degassing becomes a strongly non-linear reactive transport where volatile exsolution provides buoyancy to drive magma transport towards eruption (Stevenson and Blake, 1998;Beckett et al., 2014) while also increasing magma viscosity (Caricchi et al., 2007;Costa et al., 2009;Pistone et al., 2012) to resist flow, sometimes up to the point where shear stress and gas over-pressure lead to fragmentation and explosive eruption (see e.g. Gonnermann, 2015, and refs therein).
Experiments show that this last stage is prone to flow localization and tipping-point behaviors (e.g. Oppenheimer et al., 2015). To date, few igneous process models have included a third vapor phase. Gutiérrez and Parada (2010) included the segregation of crystals and vapour bubbles in their magma chamber model, but focus more on the thermo-chemical aspects than on mechanical coupling. Huber and Parmigiani (2018) extended the two-phase model of (Bercovici and Ricard, 2003) in the porous flow limit to two pore fluids in a compactible matrix. Afanasyev et al. (2018) employ a method based on petroleum reservoir modelling to calculate vapour and brine transport through an undeformable porous medium. However, none of these provide a model framework general enough to investigate and compare the full range of processes outlined here, including porous, mushy, and suspension flows of two, three or more solid, liquid and gas phases.

Fundamental principles 3.1 Continuum model framework
We formulate a continuum model for reactive transport in multi-phase aggregates of n material phases composed of m thermodynamic components. At each point in the continuum, we characterise phases by their system-scale velocity, v i , pressure, P i , temperature, T i , volume fraction, φ i , and component concentrations by mass, c i j . We assume that phase fractions saturate the aggregate, i φ i = 1, and component concentrations make up the entire phase, j c i j = 1. These are the n × (4 + m) primary variables in which we formulate our model.
We assume all phases are compressible, viscous materials, with densities ρ i , and viscosities, η i , that vary in space and time and are functions of independent variables. We distinguish between solids, liquids, and gases only based on their relative viscosities and densities. We approximate solids as highly viscous, nearly incompressible (∂ρ i /∂t ≈ 0) fluids. This approximation has ample precedent in geodynamics (e.g., Turcotte and Schubert, 2018), where the focus is on long-term deformation. Elastic and brittle-plastic modes of deformation are undoubtedly important for igneous processes in brittle rock. To reduce model complexity, we limit ourselves to viscous materials but point to previous work on visco-elastic/brittle-plastic two-phase flows (Keller et al., 2013;Yarushina and Podladchikov, 2015) for ways to incorporate more rheological complexity. Liquids we characterise by a relatively low viscosity and low compressibility and gases by a negligible viscosity but significant compressibility.
Our model derivation follows four steps. First, we state the fundamental principles as continuumscale conservation equations. Second, we employ thermodynamic principles to construct an expanded entropy inequality that combines all terms requiring constitutive relations. Third, we exploit the entropy inequality to choose admissible constitutive relations that describe system-scale flux and transfer processes. Whereas we choose familiar relations for intra-phase fluxes, we propose a new formulation for inter-phase transfers. Finally, we propose phenomenological closures for material response coefficients coupling fluxes and transfers to their driving forces. We thus obtain a closed model of reactive transport for n compressible, viscous phases.

Generic conservation law
We begin our model from a generic conservation law for n mass-specific vector quantities, a i , where ∂(·)/∂t is the partial derivative with time, ∇ = ∂(·)/∂x is the partial spatial derivative operator, and ⊗ is the outer vector product operator, a i ⊗ a i = a i a i T ; throughout, tensors are denoted in underlined bold face, vectors in bold face, and scalars in cursive. Unless stated otherwise, all variables and parameters are functions of spatial position and time, e.g., a i = a i (x, t). Equation (1) states that the partial phase property density, φ i ρ i a i , evolves in time due to (terms from left to right) the divergence of mass flux, q i φ ⊗ ρ i a i , carried by the volume flux, q i φ (see below for relation to phase velocity, v i ), the divergence of fluxes other than mass flux (i.e., molecular diffusion), q i a , transfers between phases, i Γ i a = 0, external sources, Q i a , and internal production, Υ i a . The production is zero for true conservation laws, but must be non-negative in the entropy equation.
Contrary to prior approaches (e.g., Drew and Passman, 2006;Bercovici and Ricard, 2003), we do not include interface stresses and surface energies in (1). The main reasons are the difficulty of accurately representing directional interface properties and hence surface stress tensors at the continuum scale (e.g., Hassanizadeh and Gray, 1990;Niessner and Hassanizadeh, 2008;Niessner et al., 2011), as well as the lack of a formalism to distinguish between motion of mass enclosed by an interface and motion of an interface with respect to the mass it encloses (e.g., Gray et al., 2013). Nonetheless, some of the leading-order effects of interface processes may later be incorporated into the phenomenological material closures.
In addition, equation (1) differs from standard continuum conservation laws (e.g., Svendsen and Hutter, 1995;Drew and Passman, 2006) in introducing the phase volume flux, q i φ . Its physical meaning and role in our formalism is best understood by revisiting the spatial averaging underpinning the derivation of continuum conservation laws (e.g., Lahey and Drew, 1988;Drew and Passman, 2006). For example, Lahey and Drew (1988) argue that the continuum-scale average of the local phase mass flux can be approximated by, where · = 1/V · dV denotes spatial averaging over the control volume, V , and the breve accent distinguishes local-scale variables from their system-scale equivalents. The approximation (2) implies that the velocity of a phase at the local scale is well represented by the average phase velocity in the control volume. That is not necessarily true, particularly in dilute suspensions. Recent theory and experiments (e.g., Segre et al., 2001;Mucha et al., 2004) show that fluctuations around the mean phase velocity can be of the same order, or even exceed the magnitude of the mean (e.g., Nicolai et al., 1995). Hence the effect of local-scale velocity fluctuations can become non-negligible (e.g., Segre et al., 2001;Mucha et al., 2004), and the volume flux carrying mass should no longer be approximated solely by the mean phase velocity, v i .
Following the reasoning of Lahey and Drew (1988) but relaxing the above assumption, we approximate the averaged phase mass flux by, The term in square brackets goes to zero only if the standard deviation of local phase velocity fluctuations, v − v , is small. Since fluctuating motion across a mean free path length leads to macroscopic diffusion, we anticipate that the phase volume flux, q i φ , in (1) will comprise a component carried on the mean phase velocity, v i , along with a second component of mechanical diffusion of phase volume as the system-scale effect of local-scale velocity fluctuations.
To write the generic conservation law (1) in Lagrangian form, we introduce the phase material derivatives in the moving reference frames of the phase volume flux, and obtain, We define the phase material derivative such that it scales with φ i . Hence, the operator vanishes if the phase is exhausted, D i φ ( · )/Dt → 0 for φ i → 0. Terms in parentheses represent the phase mass balance apart from mass transfers.

Multi-phase conservation laws
We state the conservation laws for phase and component mass, phase momentum, and phase total energy, as well as phase entropy production by substituting the conserved property densities, fluxes, transfers and sources listed in Table 1 into (1), We do not consider external sources of mass (Q i ρ = Q i j = 0). Component mass equations (6b) must sum to phase mass conservation (6a), and therefore j q i j = 0, and j Γ i j = Γ i ρ must hold. The equations (6) express the fundamental principles that underpin our model. However, they merely provide a consistent structure for enforcing conservation principles and entropy production but do not yet describe the reactive transport processes we are interested in. What is required are constitutive relations for fluxes, q i a , transfers, Γ i a , sources, Q i a , and entropy production, Υ i s , which describe the processes that may occur and the forces that are driving them.
In our constitutive choices, we adhere to three out of the four axiomatic principles discussed in Passman et al. (1984). First, the principle of local action states that thermodynamic processes respond to thermodynamic forces acting within the control volume or on its surfaces. All constitutive relations must therefore be functions only of the local state and spatial gradient of independent variables. Second, all constitutive relations must be frame invariant, that is, independent of the reference frame of an observer. The principle of frame invariance requires that, for example, processes must not depend on absolute quantities like phase velocities or their gradients, but rather on objective quantities like phase velocity differences or the symmetrical gradient tensor. And third, the entropy principle requires that constitutive relations must not violate the non-negativity of entropy production.
The fourth principle in Passman et al. (1984) states that constitutive relations for processes internal to a phase are chosen for that phase as separate from the mixture. The principle of separation of phases would require fluxes to be functions of independent variables in the respective phase only, whereas transfers could depend on variables of all phases. We do not invoke this principle because we find that processes here do not generally adhere to it. For example, if a phase is interconnected at the local-scale, then its diffusive flux can be conceptualised as separate from other phases. However, if a phase is disaggregated fluxes passing through it by necessity must intermittently pass through other phases in the control volume. Hence, fluxes must generally be considered functions of all phase variables as well.
According to these principles, we seek constitutive relations that are local, frame-invariant functions of independent variables resulting in non-negative entropy production. This admissible function space is expansive. To reduce complexity, we restrict our choice by three limiting assumptions. First, we seek only decoupled relations, where each flux or transfer is a function only of its conjugate thermodynamic force in the entropy inequality (De Groot and Mazur, 1984;Jou et al., 2001). Second, we limit our choice to linear relations, where each thermodynamic process is a linear function of its conjugate force. And third, we assume that material response coefficients coupling processes to their conjugate forces are isotropic scalars allowed to vary in space and time and to be non-linear functions of all variables. The latter assumption eliminates material anisotropy from our considerations. While there is no doubt that these limitations do not reflect the full nature of igneous systems, they represent a suitable starting point for devising models of leading-order mechanical and thermodynamic complexity.

Thermodynamic principles
We seek to assemble an expanded entropy inequality that places all transfers and fluxes explicitly under the thermodynamic constraint of non-negative entropy production (Liu, 1972;De Groot and Mazur, 1984;Jou et al., 2001). Recent works that have followed a similar procedure are the two-phase granular flow model of Monsorno et al. (2016b) , and the two-phase porous flow model of Yarushina and Podladchikov (2015). Here, we modify the approach by introducing a new, phase-symmetrical formulation for transfers between n phases. Also, we do not impose thermal equilibrium between phases (e.g., Rudge et al., 2011), but retain separate energy conservation and entropy production statements for each phase.
In analogy to the local-scale definition of total energy, we assume that specific phase total energy at the continuum scale is the sum of its internal energy, u i , and kinetic energy, 1/2 v i 2 , at the same scale: Changes in specific total energy in the phase volume flux reference frame are, Assuming that pressure-volume work is the only reversible work done on the system, u i evolves with changes to specific entropy s i , specific volume, 1/ρ i , and component concentrations, c i j : Equation (9) introduces temperatures, T i , pressures, P i , and chemical potentials, µ i j , as the thermodynamic conjugates to specific entropy, volume, and component concentrations. Allowing for other reversible work in the system would result in further contributions to (9). For example, if elastic deformation is considered terms of reversible deviatoric and volumetric strain work will be added (e.g., Yarushina and Podladchikov, 2015). Some granular flow models include a term of D∇φ i /Dt related to reversible energy in local-scale granular configurations (e.g., Monsorno et al., 2016b). These contributions present opportunities to further extend our model in the future.
Using (8) to substitute for u i , along with the equality, (9) as an expression for entropy evolution, Equation (10) now relates entropy evolution to that of total energy, momentum, and phase and component mass.
Multiplying entropy production (6e) by the absolute phase temperature, T i ≥ 0, we note that T i Υ i s ≥ 0 must hold. Hence, we write the basic entropy inequality as, We begin expanding (11) by substituting (10) for the entropy evolution term, Using conservation of phase total energy (6d), phase momentum (6c), phase mass (6a), and component mass (6b) we substitute the first four terms on the right hand-side of (12) and write, We have sorted terms according to their category as transfers, fluxes, and sources, and have used A remaining term of kinetic energy transferred by mass transfers, 1 2 v i 2 Γ i ρ , was dropped as it remains negligible in the non-turbulent flows that are the focus here.
In (13), we have interpreted the term ∂φ i /∂t as phase transfer of partial volume and hence grouped it with other transfer rates. As required of all transfers, the term sums to zero over all phases, i ∂φ i /∂t = 0. Lastly, we have introduced two terms, ±P i v i · ∇φ i , anticipating their utility in formulating frame-invariant constitutive relations. For now, we note that they have no net effect on entropy production as together they cancel out.
It is apparent from (13) that entropy in multi-phase aggregates is produced from irreversible phase interactions (transfers), from dissipative phase-internal processes (fluxes), and from interactions with the system environment (sources). The second law of thermodynamics strictly requires that entropy production of the entire system be non-negative, i Υ i s ≥ 0. One way to satisfy this constraint is to require the contribution of each of these categories to be positive in each phase (Jou et al., 2001). We note that this condition does not hold in general (see Liu, 1972) and that more rigorous treatments exist (e.g., Gray and Miller, 2014), but we nevertheless adopt it here to limit model complexity. We therefore decompose entropy production into additive contributions from transfers, (Υ i s ) trf ≥ 0, fluxes, (Υ i s ) flx ≥ 0, and sources, (Υ i s ) src ≥ 0, and subject each to the non-negativity constraint separately: We are now ready to leverage the expanded and decomposed entropy inequality (14) to choose admissible constitutive relations. Whereas our choice of fluxes will mostly follow the classical approach, we will introduce a new formulation for phase transfers. Sources will be determined by a priori knowledge of materials and the potential field environment.

Transfers
In our model, continuum-scale phase interactions comprise transfers of partial volume, mass, momentum, entropy, and energy between phases. Mass and entropy transfers act to chemically and thermally equilibrate phases. Here, we interpret momentum transfers as the phase drag resisting segregation and volume transfers as change in phase fractions due to the process commonly referred to as compaction. We use the terms segregation and compaction to refer to the deviatoric and volumetric components of phase separation. Compaction is perhaps not an optimal term since the process necessarily involves the compaction (i.e., partial volume contraction, ∂φ/∂t < 0) of some phases accommodated by the decompaction (i.e., partial volume expansion, ∂φ/∂t > 0) of others. Nevertheless, we use the term for continuity with previous work (McKenzie, 1984).
Local-scale phase interactions are driven by gradients in thermodynamic properties across interfaces. On the system-scale, these interface gradients are approximated by differences between continuum phase fields. In an n-phase model, it is expedient to express phase differences with respect to a common reference state rather than considering all possible phase pairs individually. We generically define the reference as a * = i ω i a a i , with weights i ω i a = 1 left to be determined. We assume that transfer rates are proportional to deviations, ∆a i * = a i − a * , and that they act to bring ∆a i * → 0. We emphasise that a * is not generally the state of local equilibrium but rather represents a shifting target towards which phase transfers are directed. Substituting the decomposition a i = a * + ∆a i * for the conjugate variables, v i , T i , P i , µ i j , in (14a) and collecting reference and disequilibrium terms, we write, with the added pressure-volume term decomposed as, The nonlinear term, ∆P i * ∆v i * · ∇φ i , in (16) has not previously come up in similar derivations (Bercovici and Ricard, 2003;Šrámek et al., 2007;Rudge et al., 2011). Placing it under the nonnegativity constraint implies that either ∆P i * ∆v i * must be a positive function of ∇φ i , or that it must vanish altogether. As we shall see below, the latter condition holds to a good approximation for systems with a finite phase viscosity contrast; hence, we have dropped the term. Examining terms on the first line of (15), we choose the transfer of total energy such that it does not contribute to entropy production but satisfies i Γ i e = 0, The remaining terms in (15) form conjugate pairs of transfer rates, Γ i a , multiplying their forcing disequilibria, ∆a i * . As stated above, we choose isotropic, linear, decoupled constitutive laws that satisfies non-negativity of entropy production for each conjugate pair, with non-negative, scalar transfer coefficients C i a . This form does not yet enforce the zero sum constraint, i Γ i a = 0. To do so, we choose weights, ω i a , that cast reference states, a * , as transfer coefficient-weighted sums of phase states, where we have introduced the short-hand, [( · )] Σ i , denoting normalization of an indexed quantity by the sum over index i. The reference state for each conjugate variable thus tends to the phase state associated with the dominant transfer rate coefficient. The physical intuition behind this mathematical operation is that phase transfers are driven by deviations from the state of the most rapidly equilibrating phase. Finally, our constitutive relations for transfers of entropy, component mass, momentum, and volume are, where we have normalised entropy transfers by T * and mass transfers by RT * , with the universal gas constant R. The net mass transfer is the sum of component mass transfers, Γ i ρ = j Γ i j . Additional terms ∼ ∇φ i in momentum and volume transfers ensure frame invariance as discussed in previous work (e.g., Bercovici and Ricard, 2003).
While these transfers are thermodynamically admissible and, as we will demonstrate below, bear out the known limits of two-phase flow as special cases, this is not the most rigorous or general treatment of system-scale phase interactions possible. We disregard non-linear coupling between different transfer rates and their driving forces. There is no formal reason why thermal equilibration, for example, would not generally depend on chemical potential differences, or pressure differences, etc. We also neglect the surface tensions and capillary pressures, as well as material anisotropy, all of which may be important in some regimes.
Applying transfer rates (20) to a specific problem requires closures for the rate coefficients, C i a , without which the constitutive laws remain an empty construct. We propose a possible class of closures in section 6. Further required is an equation of state for chemical potentials, µ i j . However, recognizing the substantial complexity of igneous petrology, we shall not make any attempt at it here. We note instead that it is common practice (Ribe, 1985b;Aharonov et al., 1995;Spiegelman et al., 2001;Hewitt, 2010;Liang et al., 2010;Hesse et al., 2011;Rudge et al., 2011;Keller and Katz, 2016) to simplify (20b) by introducing parameterised affinities, ∆Z i j ≈ ∆µ i * j /RT * , such that, For example, affinities can be parameterised in terms of concentration differences, Z i j = c i j −c i j , from equilibrium concentrations,c i j , prescribed as a function of pressure, temperature and composition (e.g., Keller and Katz, 2016). Alternatively, reaction rates can be eliminated from the model description by enforcing local chemical equilibrium according to parameterised phase diagrams (e.g. Katz, 2008;Weatherley and Katz, 2012;Solano et al., 2014;Jordan and Hesse, 2015;Rees Jones et al., 2018) or by employing thermodynamic phase equilibrium calculations (e.g. Tirone et al., 2009;Dufek and Bachmann, 2010;Gutiérrez and Parada, 2010;Oliveira et al., 2018). With these consideration in mind, we carry reaction rates forward by their symbol only and abstain from discussing the intricacies of igneous petrology.

Fluxes
In our framework, fluxes describe phase-internal transport of partial volume, mass, momentum, entropy, and energy. Entropy and mass fluxes act to equilibrate thermal and chemical gradients in each phase. We interpret the momentum flux as the viscous stress tensor, and the volume flux as the combination of advection on the mean phase velocity and mechanical diffusion due to local-scale velocity fluctuations.
To facilitate the choice of constitutive relations, we use the chain rule to transform vector and tensor flux divergences in (14b), with the added pressure-volume term decomposed as, We choose the flux of total energy such that it does not contribute to entropy production (terms in parentheses cancel out), The remaining terms form conjugate pairs of fluxes, q i a , multiplying their forcing gradients, ∇a i . We again choose linear, decoupled, isotropic constitutive relations that satisfy non-negative entropy production for each flux-gradient couple separately: with non-negative, scalar flux coefficients K i a . Finally, our constitutive relations for fluxes of entropy, component mass, momentum, and volume are, where we have normalized entropy fluxes by T i , and component fluxes by RT i . The use of the trace-free (deviatoric), symmetrical velocity gradient tensor, enforces symmetry of the stress tensor required to conserve angular momentum; the additional pressure and velocity terms in momentum and volume fluxes again ensure frame invariance. The component and volume fluxes constitute special cases as they are subject to constraints j q i j = 0 and i (q i φ − φ i v i ) = 0, respectively. In analogy to the technique used to enforce sum constraints on transfers, we expand forcing gradients around a reference gradient, The reference gradients are again coefficient-weighted sums, Note that generally,P * = P * , except if transfer and flux coefficients have identical relative weights, The physical intuition behind the decomposition is that diffusive fluxes subject to a sum constraint diffuse down the deviation from the forcing gradient of the fastest diffusing phase. This has the somewhat counter-intuitive effect that components/phases that are not exposed to a non-zero gradient in their phase may still experience a non-zero diffusive flux. Moreover, slow diffusers will effectively anti-diffuse as they compensate for the down-gradient diffusion of others.
As we already stated, we will not further discuss chemical potentials here. It is, however, not uncommon (e.g., Spiegelman et al., 2001;Katz, 2008;Rudge et al., 2011) to assume that ∇c i j ≈ ∇µ i j /(RT i ) and therefore to simplify chemical diffusion to, where Regarding the volume flux (26d), we note that it appears as conjugate to P i in the entropy inequality. The internally consistent choice of constitutive relation is therefore that mechanical diffusion is written as a function of phase pressure gradient deviations. To our knowledge, a constitutive relation of that form has not been proposed elsewhere. Indeed, it remains unclear how it should be reconciled with sedimentation experiments and theory (Mucha et al., 2004;Segre et al., 2001) suggesting that diffusion due to velocity fluctuations is driven by gradients in phase fractions, where and (31) implies that ∆∇φ i * ≈ ∆∇P i * /p 0 , withK i φ ≈ K i φ p 0 , and p 0 some pressure scale. We see no formal reason to rule out a diffusive flux driven by pressure gradients, nor is volume diffusion down gradients of phase fractions as well understood as thermal, chemical or momentum diffusion (see Mucha et al., 2004). For now, we retain (32) to conform with existing suspension flow models.

Sources
Source terms are determined from a priori knowledge of phase materials and the potential field environment the system is exposed to. We write the total energy source such that it does not contribute to internal entropy production (14c): For igneous systems, we take the gravitational acceleration, g, as external source of momentum, and the specific radiogenic heating rate of phase materials, H i , as external source of entropy:

Final Governing Equations
Given the constitutive relations, we now assemble the final governing equations for in the variables gives further details on how the equations are assembled, in particular on how energy conservation is transformed to a temperature equation. In the final equations, we sidestep chemical potentials by including reaction rates in symbolic form and chemical diffusion as a function of concentration gradients. We also opt for the experimentally supported form of mechanical diffusion down gradients in phase fractions.
There are (2 × n) mechanical equations, which are forms of the phase momentum and mass conservation equations arranged to emphasise multi-phase interactions: Here, we have introduced a short-hand for phase velocity differences or "segregation velocities", v i ∆ = φ i ∆v i * , and phase pressure differences or "compaction pressures", P i ∆ = φ i ∆P i * . These generalise the Darcy's segregation velocity (Hubbert, 1957) and Terzaghi's effective pressure (Skempton, 1960), both commonly used in porous flow models. Given the constitutive choices above, equations (36) express that segregation is driven by velocity changes (acceleration), gradients in reference and compaction pressures, momentum diffusion (shear stresses), momentum transfer by reaction, and buoyancy, and is modulated by the inverse momentum transfer coefficient, Compaction is driven by density changes (compressibility), divergence of reference and segregation velocities, mechanical diffusion, and volume transfer by reaction, and is modulated by the inverse volume transfer coefficient, Summing (36) over all phases yields the mixture momentum and mass conservation, written for the gradient of the mixture pressure,P = P * + i φ i ∆P i * , and the divergence of Alongside the mechanical equations, we write thermo-chemical evolution as (n) equations for phase fractions, (m × n) for component concentrations, and (n) for phase temperatures, where α i is the phase thermal expansivity, c i P the phase specific heat capacity, and L i the latent heat of phase change. According to our constitutive choices, equations (38) express that phase fractions and component concentrations evolve due to advection, diffusion and reaction. The evolution of temperature (sensible heat) is driven by thermal diffusion, thermal equilibration, latent heat of reactions, adiabatic heating, radiogenic heating, and the dissipation of deviatoric, segregation, and compaction deformation.
These final governing equations present a physically consistent and thermodynamically admissible model description of multi-phase reactive transport at the continuum scale. However, without closures for flux and transfer coefficients and equations of state for density and other material properties, these equations remain an empty structure. We will not discuss equations of state for igneous materials here, but will proceed to consider the problem-specific coefficient closures required to fill this mathematical construct with physical meaning.
6 Phenomenological Closures 6.1 Diffusive transport In the constitutive relations above, fluxes and transfers are related to their conjugate forces by material response coefficients, which we have limited to non-negative, scalar functions of system variables. To guide our choice of coefficients, we can no longer rely on fundamental principles. Instead, it is necessary to invoke a specific natural or engineering problem, and to rely on empirical evidence and conceptual thought in choosing non-unique, phenomenological closures for coefficients modulating fluxes, K i a , and transfers, C i a . The class of closures we propose here is based on the general assumption that transfers and fluxes are facilitated by diffusive transport within and between phases. Our material closures will therefore depend on the diffusive properties of the pure phase materials. In addition, they must contain information regarding the system-scale effect of local-scale phase distribution and topology in the aggregate.
Regarding the pure-phase diffusive properties, the thermal diffusivity, κ i T , chemical diffusivity, κ i j , and momentum diffusivity, κ i v = η i /ρ i (with κ i v the kinematic, η i the dynamic viscosity) of igneous materials are relatively well known from experiments and theory. The mechanical diffusivity, κ i φ , is a concept developed based on suspension flow experiments and theory (e.g. Segre et al., 2001;Mucha et al., 2004). As all diffusivities, it is understood as the product of a mean velocity deviation times a mean free path length (or correlation length), Whereas possible models for mean local-scale phase velocity fluctuations, δv i , have been discussed in significant detail (see Tee et al., 2002;Mucha et al., 2004), we will limit ourselves to the leadingorder approximation, δv i ∼ |∆v i * |. The correlation length is thought to be 2 − 20 × the suspended particle size (e.g., Segre et al., 2001). In the dilute suspension limit, |∆v i * | between the suspended and carrier phases is well approximated by the magnitude of the hindered Stokes settling speed, with ∆ρ i * = ρ i − k ω k v ρ k the density contrast to the carrier phase and a hindering exponent of m = [4, 6]. In the pure-phase limit, we therefore approximate the mechanical diffusivity as, To give an example, in a dilute suspension of 3 mm sized crystals with a density constrast of 500 kg/m 3 to a carrier melt of 10 3 Pas viscosity, the Stokes settling speed is ∼ 5 × 10 −5 m/s. Assuming a correlation length of ten times the granular size, the resulting mechanical diffusivity is κ i φ ≈ 10 −6 m 2 /s, of same order as a typical thermal diffusivity.
Regarding the system-scale effects of local-scale phase topology, we assume that they may, to a leading order, be captured by two phenomenological attributes. The first is the characteristic length scale, d, representing the mean size of phase constituents (e.g., mean grain size). The second is a metric termed transmission function, X ik φ , we introduce to describe the transition between connected and disconnected phase topologies. Transmission functions for each phase have n components: Intra-phase transmissions, X ik φ , i = k, describe the availability of connected pathways along which diffusive transport is transmitted; inter-phase transmissions, X ik φ , i = k, express contact with other phases through which diffusive transport is transmitted to adjacent phases surrounding isolated constituents of a disaggregated phase. These transmission functions are conceptually related but not identical to geometric attributes such as phase contiguity, relative contact area, volumetric connectivity, etc., which have been theoretically derived for idealised equilibrium melt-grain textures (e.g., Takei, 1998;Takei and Holtzman, 2009;Zhu and Hirth, 2003;Ghanbarzadeh et al., 2014;Rudge, 2018a), or empirically quantified by microtomographic analysis of analogue systems (Zhu et al., 2011;Miller et al., 2014;Colombier et al., 2018). Generalizing on this range of phase topology metrics, phase connectivities, x ik φ , are generally defined as fractions of total surface area of phase constituents of phase i touching other constituents of the same phase, x ii φ , or another phase k, x ik φ . Our transmission functions relate to these phase connectivities as X ik φ ∼ 1 + log 10 x ik φ . Generally, phase connectivities, and hence the transmission of diffusive transport within and between phases, may depend on factors such as phase fractions, topology, size distribution, and relative surface energies, and the rate and state of deformation. Here, we intentionally bypass the potential complexity of modelling phase connectivities and instead directly prescribe phenomenological transmission functions of the form, The three fitting parameters are, slopes A i that describe a gradual decrease in phase connectivity with decreasing phase fractions; critical phase fractions 0 ≤ B ik ≤ 1 ( k B ik = 1) that specify disaggregation or percolation thresholds towards which phase connectivity drops off steeply; weights 0 ≤ C i ≤ 1 set the relative smoothness of that step function that models that phase disconnection. Figure 3 shows an example of transmission functions for an igneous two-phase, solid-liquid system ranging from partially molten rock and partially crystalline magma. The solid intra-phase transmission, X ss φ , shows where a contiguous granular matrix is present, while that for the liquid, X φ , marks where melt forms a connected phase. The liquid may form a connected network along grain boundaries even where the solid is contiguous. The inter-phase transmissions, X s φ and X s φ , grow towards unity where the solid and liquid phase become disaggregated grains and disconnected melt pockets, respectively. The reference fitting parameters are chosen for a solid disaggregation threshold at 0.30, and a liquid percolation limit at 0.02 melt content. Parameter values for the reference curves are given in Table A1.
We have calibrated the reference example transmissions in Fig. 3 to reflect a scenario of a low dihedral angle (∼ 30 • ) between melt and rock phases, as for basaltic melt in an olivine-rich rock (e.g., von Bargen and Waff, 1986). For comparison, we show the transmission-equivalent of granular contiguity (black line, 1 + log 10 ϕ s /16) calculated from the semi-empirical formula, ϕ s = 1 − 2φ 0.5 (Takei and Holtzman, 2009). Our solid transmission not only captures the loss of contiguity as the solid disaggregates (contiguity drops to zero), but additionally describes a smooth transition from intra-phase connectivity within the solid to inter-phase connectivity between solid and liquid.
Since in our framework, transmissions prescribe the effective multi-phase behavior of a specific natural or engineering system, their calibrating becomes the most consequential step of applicationoriented model building. While the calibration remains phenomenological and non-unique, it presents a flexible tool for transparently testing specific hypotheses regarding these complex natural systems.

Flux coefficients
In our framework, fluxes capture the system-scale effects of local-scale diffusive transport within each phase. We seek to formulate phase-wise flux coefficients that express the degree to which effective transport is permitted based on the properties of a phase and the leading-order attributes of its local-scale topology. We do so by introducing flux permission functions, θ i a ≥ 0, for each phase that are dimensionless, non-negative functions of pure-phase diffusivities and the phase transmission functions introduced above. Permissions are equal to unity for fluxes passing along straight, connected pathways of a phase and increase or decrease across a disaggregation or percolation threshold depending on whether disconnected phase constituents are now connected to other phases of higher or lower diffusivity.
To obtain the correct physical dimensions according to (6), we now choose flux coefficients of the form, Thermal and chemical diffusivities in igneous materials do not tend to vary as strongly as phase viscosities, on which both our momentum and volume flux coefficients depend. Hence, the following considerations are most pertinent to momentum and volume fluxes, which depend the effective rheology of a disaggregating solid (Arzi, 1978;Van der Molen and Paterson, 1979;Renner et al., 2000;Costa et al., 2009), or a liquid disconnecting at a percolation threshold (von Bargen and Waff, 1986;Zhu and Hirth, 2003).
The study of the effective rheology of aggregates has a long tradition in Effective Medium Theory (e.g., Mavko et al., 2009). Effective Medium Theory has mostly been concerned with elastic properties (Hashin and Shtrikman, 1963;Hill, 1963;Watt et al., 1976); it casts effective deformational properties as weighted averages of pure-phase properties with weights that take into account relative proportions and topologies of local-scale phase constituents, and averaging techniques that reflect the partitioning of flux and forcing gradients between phases. For example, if the total flux accommodated by the aggregate is the sum of partial phase fluxes, while all phases experience the same forcing gradient, a weighted arithmetic average is used. This provides an upper limit to the effective flux coefficient known as the Voigt bound (Hill, 1963). Or, a harmonic mean is called for if the total forcing gradient supported by the aggregate is equal to the sum of partial phase gradients, while each phase sustains the same flux. This lower limit is known as Reuss bound (Hill, 1963). Stricter bounds have been derived by Hashin and Shtrikman (1963) for an idealised scenario of spherical inclusions. Yet, these bounds still leave considerable room for different effective behavior in between.
For a given set of weights, geometric averaging produces effective properties halfway between upper and lower bounds in logarithmic space, a convenient property when dealing with possibly tens of orders of magnitude phase viscosity contrasts. We therefore write permission functions as weighted geometric averages of pure-phase diffusive flux coefficient contrasts, where the weights are power laws of the transmission functions, The two dimensionless fitting parameters, a prefactor D ik and an exponent E ik , allow flexible scaling of permission functions to problem-specific contexts, for example to fit empirical data or canonical models in specific limits such as porous or suspension flows. Figure 4: Sample family of momentum and volume flux permission and effective flux coefficient as a function of liquid fraction for example igneous two-phase system with solid and liquid viscosities of 10 1 8 Pas and 10 2 Pas, respectively, granular length scale of 3 mm and phase density contrast of 500 kg/m 3 . All curves are based on the reference transmission curves in Fig. 3. The shape of momentum flux permissions (a) and volume flux permissions (b) for the solid (blue) and liquid (red) phases reflect that of the underlying transmission curves, whereas the amplitude is controlled by phase viscosity contrasts. Resulting coefficients for momentum (c) and volume (d) fluxes show the variability prescribed by permissions. Comparison of the reference coefficients (heavy) with previous models of effective viscosity (e) and mechanical diffusivity (f ) shows good correspondence.
Returning to our example two-phase system above (Fig. 3), Figure 4 shows the momentum and volume flux permissions (see Fig. 4(a) & (b)) and the corresponding flux coefficients (see Fig. 4(c) & (d)). All curves in Fig. 4 are based on the reference phase transmissions in Fig. 3 and on a granular length scale of d = 3 mm, pure-phase rock viscosity of η s = 10 18 Pas and melt viscosity η = 10 2 Pas, and a phase density contrast of ∆ρ sl = 500 kg/m 3 . Fitting parameters for the reference curves are given in Table A1. Owing to the strong phase viscosity contrast, permissions vary over sixteen orders of magnitude, and with much of that variation occurring over a limited region of two-phase space.
The solid momentum flux permission ( Fig. 4(a)) describes the variation of effective solid rheology. After an initially log-linear decrease it steeply drops across the disaggregation limit. With the liquid transmission increasing sharply from the percolation threshold, the liquid momentum permission first steeply plunges from a high value that reflects the effective rheology of isolated melt pockets before gradually approaching unity in the pure-melt limit. The permissions and flux coefficients for mechanical volume diffusion show inverse trends to momentum diffusion since phase contrasts in mechanical diffusivity reduce to the inverse ratio of phase viscosities, M ki The example illustrates key aspects of effective momentum and volume fluxes transitioning between the end-members of porous and suspension flows in igneous two-phase systems. Figure  4(e)-(f) shows a comparison of our flux coefficients with previous work. As expected, the viscous equivalents of the Hashin-Shtrikman bounds describe a wide envelope around the effective solid and liquid viscosities produced by our model (η i eff = η i θ i v ). To benchmark the increase in liquid viscosity by the addition of rigid particles we use the Einstein-Roscoe relation, η eff = η (1 − φ /φ c ) (φcB) (Einstein, 1906(Einstein, , 1911Roscoe, 1952;Krieger and Dougherty, 1959). Our example produces an effective liquid viscosity corresponding to an Einstein-Roscoe curve with φ c = 0.58 and B = 2.5.
On the solid-dominated end of two-phase space, a number of effective shear viscosity laws for the disaggregating solid have been proposed. A widely adopted empirical relation describes an exponential decrease with melt fraction, η s eff = η s exp(−λφ ) (e.g. Hirth and Kohlstedt, 2003). We compare our effective solid viscosity to an exponential law with slope λ = 27 (Mei et al., 2002). Furthermore, we compare our fit to the theoretically derived relations of Takei and Holtzman (2009) and Rudge (2018b), who based their closures on calculations of idealised grain-melt equilibrium textures. Our example first follows a similar curve to that of Rudge (2018b), but continues on to a smooth step across the disaggregation threshold rather than dropping to zero. The curve of Takei and Holtzman (2009) shows a similar initial slope and disaggregation drop-off apart from the slight offset at the onset of melting, a peculiarity of the grain boundary-diffusion model they applied. The fitting model of Schmeling et al. (2012) uses relations from effective media theory to produce similar curves.
In order to connect effective rheology models between the porous and suspension limits, Costa et al. (2009) proposed a semi-empirical fitting law to calibrate effective aggregate viscosity curves across the entire two-phase space. Since their formalism has the same general trends in mind as ours, we find it straightforward to calibrate one of their curves to closely correspondes to our reference example (parameter values listed in Table A1).
Comparing our mechanical diffusivity closures to previous work, we first return to the approximate law relating velocity fluctuation magnitudes in the suspension limit to the hindered Stokes settling speed eq:hindered-Stokes. Our closure for the solid mechanical diffusivity gives a reasonably good fit given a correlation length h s = 10d s , and a hindered Stokes settling speed with m = 5. We further find that the theoretically empirically confirmed relation of Segre et al. (2001) taking into account effective liquid viscosity and structure factors derived from Brownian suspension theory (e.g., Al-Naafa and Selim, 1993) produces a slightly steeper but still very similar curve.

Transfer coefficients
In our framework, transfer rates capture the system-scale effects of local-scale phase interactions. We include entropy transfers governing thermal equilibration, component mass transfers governing chemical equilibration, momentum transfers governing segregation, and volume transfers governing compaction. We assume that transfers are limited by inter-phase diffusive transport over the local length scale, d. Note that advective transfers in our framework are captured as properties carried by net mass transfers. For example, v i Γ i ρ is the advective transfer of momentum. However, Γ i ρ is the sum of component transfers, Γ i j , which we again assume to be limited by chemical diffusion across phase constituents. Hence, we write closures for entropy, component, momentum and volume where we have used the granular-scale buoyancy pressure p 0 = ∆ρ i * gd to obtain the volume transfer coefficient, C i φ , in units of (Pas) −1 . We scale transfer coefficients with φ i (1 − φ i ) to ensure that transfers to and from a phase must cease both in the pure-phase (φ i → 1) as well as the phaseexhaustion (φ i → 0) limits.
As for fluxes, we have introduced dimensionless, positive functions ψ i a as transfer permissions to capture the variation of effective transfer rates as a function of phase diffusivity contrasts and phase transmissions, with a fitting prefactor, F ik , and an exponent, G ik , allowing flexible scaling to fit problem-specific behavior.
We will again focus the following considerations around momentum and volume transfer permissions. Owing to strong phase viscosity contrasts, they are expected to show more variability than thermal and chemical ones. The two related processes of segregation and compaction also make up the most consequential phase interactions controlling dynamics in our framework.
We define segregation as the relative motion of one phase in the aggregate with respect to the others: the deviatoric component of phase separation. Momentum transfers govern segregation: fast transfer tends to dampen segregation, whereas slow transfer tends to permit velocity differences to persist or grow. Segregation mobility therefore is inversely related to the rate of momentum transfer. According to the canonical models, the segregation mobility of a pore liquid is limited by its own phase viscosity, whereas the mobility of suspended phase is limited not by its own viscosity but by that of the carrier melt. We will now show how our proposed momentum transfer coefficients reconcile these contrasting behaviors and bridge the two limits of porous and suspension flows.
Returning to our two-phase example, we find that both end-member behaviors are captured well by momentum transfer permissions very similar in shape to the flux permissions above, ψ i v ∼ θ i v . Figure 5 shows momentum transfer permissions (a) along with their resulting transfer coefficients (c). These are calculated from the same reference transmission functions and phase properties as flux permissions and coefficients in Fig. 4. Fitting parameters used for the reference curves are listed in Table A1. Notably, the permission exponents, G i , chosen to produce the desired behavior, are the same as the ones chosen for flux permissions, E i , above. The flux and transfer permission prefactors, D i and F i , are within about an order of magnitude.
To put momentum transfer coefficients into context with previous work, we compare momentum transfer coefficients to the canonical Darcy coefficient: where k i φ is the Darcian permeability. In Fig. 5(e), we plot the segregation coefficients (right hand side of (48)), corresponding to the reference momentum transfer coefficients in Fig. 5(c). We (1 − φ ) −2 , in the porous flow limit, as well as a bracket of hindered Stokes coefficients, in the suspension flow limit. In both limits, our approach produces segregation coefficients that fit well with expected behavior without resorting to any additional model assumptions in addition to the chosen transmission functions and pure-phase material properties.
Our approach provides a new perspective on compaction. We define compaction as the evolution of phase fractions in response to phase pressure differences: the volumetric component of phase separation. In his seminal work on porous melt transport, McKenzie (1984) introduced a volumetric stress term, ζ s ∇· v s , to the matrix shear stress tensor. The material parameter, ζ i , denoted "bulk viscosity", has since been widely discussed in that context (e.g., Simpson et al., 2010;Schmeling et al., 2012;Rudge, 2018a). However, the term may be considered misleading as it suggests a parallel to either effective aggregate properties such as the bulk density or else the bulk modulus in elastic theory. In the latter sense, the bulk viscosity would describe a phase material's resistance to irreversible volumetric deformation in response to pressure, ∇·v i = ζ i −1 P i . In contrast, compaction is fundamentally a multi-phase process and relates transfers of partial volume between phases to phase pressure differences, ∂φ i /∂t ∼ ∆P i * (e.g., Bercovici and Ricard, 2003).
Existing porous flow models mostly formulate compaction in terms of the solid-liquid pressure difference, φ (P s − P ), also known as Terzaghi's effective pressure. Though originally introduced in the context of soil failure (Terzaghi, 1943;Skempton, 1960), the same definition has also been employed as compaction pressure in models of porous flow in water-saturated sediments (Birchwood and Turcotte, 2012;Morency et al., 2007) and partially molten rock (Connolly and Podladchikov, 1998;Katz et al., 2006;Keller et al., 2013). Whereas McKenzie's compaction model implies that ∇ · v s is the only process contributing to ∂φ /∂t, and that the solid matrix is the load-bearing phase, Bercovici and Ricard (2003) extended the two-phase model to a phase-symmetrical form, andŠrámek et al. (2007) further included reactions as an additional contributor to compaction rate. Here, we expand on these models and present an n-phase, phase-symmetrical formulation, which relates generalised phase pressure differences, ∆P i * , to volume transfers, ∂φ i /∂t, that comprise compressibility effects and mechanical volume diffusion in addition to divergent phase flow and reactions.
To minimizing the number of separate model choices, we exploit model symmetry and propose volume transfer coefficients (46d) depending on the same mechanical diffusivity (41) as volume flux coefficients above. While this formulation seems to imply a physical link between phase volume transfers and fluctuation-driven volume diffusion, we do not see any reason to suspect such a link exists beyound the shared dimensionality between the terms. Still, we find that by scaling κ i φ by the granular-scale buoyancy pressure, p 0 = ∆ρ i * gd 2 , while applying the same formulation as for other transfer coefficients (C i a ∼ κ i a /d 2 ), the result of C i φ ∼ φ i (1 − φ i )/η i recovers a form equivalent to the compaction viscosity employed in previous work (McKenzie, 1984;Bercovici and Ricard, 2003;Simpson et al., 2010), with the two being inversely related, ζ i = φ i 2 /C i φ . We use this latter form to define compaction coefficients relating to volume transfer coefficients analogously to how segregation coefficients relate to momentum transfer coefficients in (48).
Suspension flow models typically do not include a compaction term. Rather, they implicitly assume that pressure differences vanish in that limit. An exception to the rule is Bürger et al. (2000), who derive a coupled sedimentation-consolidation model, but employ a different approach to compaction we do not find straightforward to reconcile with other work. Even though we find vanishing compaction pressure in the suspension limit a reasonable approximation, we wish to impose minimal assumptions and will therefore maintain model symmetry by retaining the volume transfer terms for all phases in all limits.
With these considerations in mind, we now return to our two-phase example. Figure 5 shows volume transfer permissions (b) along with their corresponding volume transfer coefficients (d). Again, we note that volume transfer permissions follow inverse trends compared to momentum transfer permissions, owing to M ki φ ∼ M ki v −1 . Fig. 5(f) shows that the corresponding compaction coefficients follow the general trend of the effective shear viscosity curves in Fig. 4(e). In the high-melt fraction limit, the liquid compaction coefficient exceeds that of the solid phase, taking on the role of the load-bearing phase for compaction from the solid phase. Our closures result in compaction viscosities tending to infinity in both porous and suspension limits, thus fulfilling the important role of enforcing incompressible flow for incompressible phase materials in the pure-phase limits.
In previous work, the solid compaction viscosity is generally assumed to relate to shear viscos-ity as ζ s ∼ η s /φ (McKenzie, 1984;Bercovici and Ricard, 2003;Simpson et al., 2010;Schmeling et al., 2012). In contrast, Takei and Holtzman (2009) find no inverse dependence on φ for grain boundary-diffusion creep in idealised grain-melt textures, and Rudge (2018b) obtains a less pronounced dependence of ∼ 1/ log φ for the case of vacancy-diffusion creep. Our closures take noticeably higher values at low melt fractions than these two theoretical studies suggest. How the results of Takei and Holtzman (2009) and Rudge (2018b) would generalise to non-idealised grain geometries or to other creep mechanisms remains unclear. However, comparing our coefficients to the exponentially weakened shear viscosity of Hirth and Kohlstedt (2003) divided by φ , and the fitted aggregate viscosity after Costa et al. (2009) divided by φ (1 − φ ) shows that our closures recover the canonical inverse dependence on phase fractions. In any case, our formulation is built such that it may be calibrated to a wide range of compaction-to-shear viscosity ratios.

Discussion
With the governing equations and material closures in place, we can now use our model framework to construct fully closed model cases and to discuss some important ramifications for understanding reactive transports in igneous multi-phase systems. The broad range of flow regimes from the porous to the suspension limit are difficult to capture in current model frameworks, particularly when exsolution of magmatic volatiles adds a third phase. Our model provides a more general modelling vantage point from where to comparatively study the diversity of reactive transport environments arising in igneous systems from source to surface. In the following, we will first demonstrate how to recover the known special cases of two-phase porous and suspension flows of solid rock or crystals, (i = s), and liquid melt (i = ). We will then add a third phase of compressible volatile phase (i = v), where the latter may be a supercritical fluid, a magmatic brine, or a gaseous vapour. For simplicity, we will do so under the assumptions of local thermal equilibrium (T i = T ) and an isochemical aggregate (c i j = c 0 ). To include reactions without specifying a thermodynamic model for chemical potentials, we assume a known melting rate (mass transfer from from solid to liquid phase), Γ s ρ , and degassing rate (mass transfer from liquid to vapour phase), Γ v ρ .

Two-phase porous to suspension flows
The governing equations for low-Reynolds number, incompressible, reactive two-phase flows are an important special case of the general equations (36)- (38): The two momentum equations, (49a) & (49b), are equivalent to the phase-symmetrical, two-phase momentum equations given by Bercovici and Ricard (2003) when neglecting surface tensions. The two-phase drag coefficient, c, in Bercovici and Ricard (2003) relates to our momentum transfer coefficients as, The two mass equations, (49c) & (49d), differ formally from Bercovici and Ricard (2003) in that they are not written as evolution equations for density or phase fractions, but rather as statements for phase compaction pressures. Our formulation also includes the mechanical diffusion terms their theory (Bercovici et al., 2001;Bercovici and Ricard, 2003) did not include. In analogy to the twophase drag coefficient, our model suggests a corresponding two-phase volume transfer coefficient, Whereas we have retained the phase evolution equation for the solid phase (49e), the solid and liquid phase equations are exactly equivalent in saturated two-phase systems where φ s = 1 − φ . In the temperature equation (49f), we have assumed constant and equal heat capacity, c P , thermal expansivity, α 0 , and internal heating rate, H 0 between the phases.K s = i K i s is the aggregate entropy flux coefficient (thermal conductivity),ρ = i φ i ρ i the aggregate density, L is the latent heat of melting, andΨ d is the total viscous dissipation of shear, segregation and compaction flow. The aggregate material derivative in the temperature equation written in terms of reference and segregation velocities becomes, Lastly, we have assumed that pressure changes in the adiabatic term in (49f) are dominated by motion against gravity, To write reduced equations in the limits of small and large melt fractions, we consider the physical interpretation of the reference velocity and pressure, v * and P * , and consequently the segregation velocities and compaction pressures, v i ∆ and P i ∆ . Returning to our two-phase example from above, Fig. 6 shows the reference velocity and pressure weights calculated from the relative weight of momentum and volume transfer coefficients. Figure 6: Porous to suspension flow transition for example two-phase system. Reference weights for velocity (a) and pressure (b) for the solid (blue) and liquid (red) as function of liquid fraction. Segregationcompaction lengths for liquid-in-solid (red) and solid-in-liquid (blue) shown in (c). Dotted lines indicate liquid percolation (φ = 0.02) and solid disaggregation (φ = 0.3) thresholds. Porous (blue), Mush (yellow), and Suspension (red) regimes indicated as shaded regions.
In the porous flow limit, the solid velocity weight, ω s v , asymptotically approaches unity, whereas ω v vanishes. As a result, v * → v s for φ → 0, meaning that the solid velocity becomes the reference point for segregation. In contrast, we find the opposite trend for pressure weights, ω i φ . Hence, P * → P , meaning that the liquid pressure becomes the reference point for compaction. Both sets of weights are a consequence of the momentum and volume transfer permissions being a function of phase viscosity contrasts and their inverse, respectively. Based on conceptual reasoning, Bercovici and Ricard (2003) anticipated similar weights for what they termed interfacial velocity and pressure. Their equivalent weights, , are indeed the same as ours if transfer permissions are identical to unity. Their pressure weights were simply given as ω i φ = 1 − ω i v , which leads to the same result as ours in the porous limit. Whereas we introduced our weights to impose the zero sum constraint on transfer rates above, we find that the result introduces physically relevant reference states for transfer processes in an n-phase context.
Examining the phase segregation velocities and compaction pressures in the porous limit, we find that v s ∆ and P ∆ vanish. Our choice of transfer coefficients therefore implies that the solid velocity, v s , the liquid pressure, P , the liquid segregation velocity, v ∆ , and the solid compaction pressure, P s ∆ are suitable variables to pose the reduced porous flow equations. The latter two variables are of course identical to the Darcy segregation velocity, v ∆ = φ (v − v s ) (Darcy, 1856;Hubbert, 1957), and the Terzaghi effective (or compaction) pressure, P s ∆ = φ s (P s − P ) (Terzaghi, 1943;Skempton, 1960). Whereas previous models have cited empirical thought or tradition to justify these particular variables to write their porous flow models, our framework formalises that choice as a direct consequence of chosen phenomenological closures for phase transfer coefficients. With this, write the governing equations in the porous flow limit as, ∇· In (53a),K v = i K i v is the aggregate momentum flux coefficient (i.e., effective aggregate viscosity). The way we have arranged the reduced two-phase equations highlights that porous flow in a deformable matrix is governed by Stokes flow of the mixture (53a)-(53b), coupled to the compaction of the solid matrix (53d) and the segregation of the pore liquid (53c). Considering the relatively small magnitude of liquid momentum and both solid and liquid volume flux coefficients in this limit, (see Fig. 4(e), (f)), we have dropped shear stress terms for the liquid as well as mechanical volume diffusion terms for both phases. The aggregate momentum and mass equations, (53a)-(53b) recover the incompressible Stokes equations if the segregation velocity, compaction pressure and melting rate vanish. The third and fourth equations recover the explicit relations for segregation velocity and compaction pressure according to Darcy's law for segregation and McKenzie's law for compaction, but now written in terms of momentum and volume transfer coefficients. Substitution of (53c) and (53d) into the right-hand sides of the aggregate equations (53a)-(53b) would yield the original form of McKenzie (1984).
For phase evolution (53e), we retain the solid equation, in which the divergence of v s ∆ drops out, leaving only reference velocity terms. In the temperature evolution (53f), we have spelled out non-negligible contributions to viscous dissipation, which are those of solid shear, liquid segregation, and solid compaction flow. Equation (53b) marks a slight but potentially consequential shift from previous models of reactive porous flow (Katz, 2008;Rudge et al., 2011;Keller and Katz, 2016;Turner et al., 2017). In our model, the compaction pressure depends on both the compaction rate (∇· v * ) and the volume change of reaction (Γ s ρ /ρ s ), whereas previous models had related it to the compaction rate only (e.g., Katz, 2008;Keller and Katz, 2016). Hence, at steady-state (D * φ i /Dt = 0), compaction and reaction rates in our model must balance out, and the compaction pressure vanishes. This may have ramifications for melt focusing beneath mid-ocean ridges. Turner et al. (2017) recently proposed that melting/crystallization reactions produce a compaction pressure component, P s ∆ ∼ ζΓ/ρ s , which remains at steady state. The resulting melting rate-induced pressure gradients were found to be potentially greater than those due to phase buoyancy contrast. As a result, Turner et al. (2017) reported that this mechanism significantly aided melt focusing from the ridge flanks towards the center of the melting regime. Our model suggests that the mechanism may need to be reevaluated.
Scaling analysis of the porous flow equations reveals three velocity scales that capture the zeroth-order flow dynamics: the aggregate convection speed, u * conv , liquid segregation speed, u segr , and solid compaction speed, u s comp . To write these, we scale pressures as [P * , P s ∆ ] ∼ ∆ρ 0 g 0 0 , with ∆ρ s 0 the phase density difference scale, g 0 the gravity, and 0 a length scale characteristic of the system size (e.g. domain depth). We further scale velocities as WithK v,0 ∼ η s 0 θ s v,0 , the convective speed becomes the Stokes speed of a viscous diapir of radius 0 , driven by the buoyancy difference, ∆ρ 0 g 0 , and resisted by the solid viscosity, η s 0 . The two transfer coefficient scales expressed in terms of viscosity and permission scales are C v,0 = φ 0 η 0 ψ v,0 /d 2 0 , and C s φ,0 = ψ s φ,0 /η s 0 , with the granular length scale, d 0 . All permission scales, θ s v,0 , ψ v,0 , ψ s φ,0 are taken as their functional value evaluated at the scaling phase fraction, φ 0 . Two dimensionless numbers useful to characterise flow regimes are found by taking ratios of the segregation and compaction speed scales by u * conv : The first we term the segregation number, R segr , which is a function of the effective phase viscosity ratio and the square ratio of the granular by the system length scales. The second, which we term compaction number, R s comp , reduces to the product of the solid momentum flux and volume transfer permission scales; it corresponds to the solid shear-to-compaction viscosity ratio, η s /ζ s , in canonical notation. The two numbers can be combined to find the segregation-compaction number,

Continuum model of igneous reactive transports
This dimensionless number can be interpreted as the squared ratio of a segregation-compaction length scale, δ s s/c , by the system length scale, 0 . The former is the inherent physical length scale associated with the liquid phase segregating through the compacting solid phase, and has been identified and termed the compaction length in previous work (McKenzie, 1984, and following works). Here, we generalise the scale to the broader multi-phase context. For n-phase systems, several such length scales may arise, each being the length scale characteristic of one phase i segregating through another phase k, which compacts to accommodate the relative motion of the first, Consistent with this definition, we generalise the segregation and compaction numbers for any phase i with respect to convection of the aggregate, The segregation-compaction length ("s/c-length" for short) of phase pairs is an effective tool for assessing the leading order multi-phase flow behavior. Since R i comp is typically close to unity or larger, the s/c-length can be used as a combined diagnostic for the length scale of segregation as well as its relative importance over both compaction and convection in a system of a known size. In Fig.  6(c), we show δ s s/c (liquid segregating through compacting solid, red line) as well as its counterpart, δ s s/c (solid segregating through compacting liquid, blue line), both truncated below 1 mm. At low melt fractions, the liquid-in-solid s/c-length rapidly increases from the granular scale to order ten meters just above the percolation threshold at φ = 0.02. We define this steep initial increase as the onset of the porous flow regime (blue shading in 6). The s/c-length continues to increase until it eventually saturates above order 1 km at φ = 0.12. Up to this point, the addition of further melt increases the speed of segregation relative to both convection of the aggregate and compaction of the solid, and also increases the length of connected pore space communicating compaction pressure gradients through the aggregate. Qualitatively, therefore, melt extraction grows more efficient and melt drainage systems grow more extensive as more melt is added in the porous flow regime. At intermediate melt fractions, δ s s/c begins to contract again as the disaggregating solid phase looses strength. Addition of melt now leads to solid phase creep picking up relative to melt segregation. Also, the contraction of connected pore space with increasing melt content implies that heterogeneities prompting melt accumulation may now incite the formation of localised melt-rich lenses. If melt has lower density than the surrounding mush body, such lenses may become gravitationally unstable, creating melt-rich diapirs rising through crystal-rich (Cashman et al., 2017;Seropian et al., 2018). We define the turning point from growing to contracting segregation/compaction length as transition from the porous to the mush regime (yellow shading in 6). In our example, that transition occurs at φ = 0.12. From there, δ s s/c drops gradually back down to the granular scale and indeed falls below the value of its counterpart, the solid-through-liquid s/c-length, δ s s/c . We choose the cross-over point occurring at φ = 0.72 in our example to mark the transition from the mush to the suspension flow regime (red shading in 6).
Our model framework suggests that the mush regime is particularly challenging to model. It requires solving the phase-symmetric equations (49) for all phase velocities and pressures since the liquid and solid phase now contribute similarly to the dynamics, which we see in the velocity and pressure weights taking on intermediate values. A further complication is the drastically varying inherent length scale from several km down to the grain scale over an increase of less than 20% in melt fraction. Our analysis suggests that local-scale processes including flow localization and interface formation will become increasingly important with growing melt fraction in the mush regime, which are challenging to incorporate into a mixture-theory framework such as the one adopted here. Clearly, more work is needed to fully characterise dynamics in this important transitional regime.
In the suspension flow limit, the liquid velocity weight increases to dominate the reference velocity, although the partitioning is not as clear as in the porous regime. Still, it holds to within a reasonable approximation that v * → v , and v ∆ ≈ 0. Hence, we formulate the reduced equations for two-phase suspension flow in the velocity variables v * and v s ∆ . As regards pressure, the reference weights of both phases approach 0.5 at φ s → 0, and thus pressures equally contribute to dynamics. However, with the inherent length scale dropping to just above the granular scale, phase pressure differences scale as ∼ ∆ρ 0 g 0 d 0 , predicting negligible pressure differences of 1 Pa for typical material parameters. Hence, to a good approximation, pressure differences vanish, and the governing equations may safely be written for a single pressure, P * = P s = P . With this, write the governing equations in the suspension flow limit as, Here, we have expressed all phase fractions in terms of the solid fraction, φ s , and used the phase evolution equation for the liquid phase to write (59d). We assumeK φ ≈ 0 in this limit consistent with v ∆ ≈ 0, and thus the mechanical diffusion flux depends on solid mechanical diffusivity and gradients in solid fraction only. Our suspension equations broadly agree with previous theory (e.g., Drew, 1971;Monsorno et al., 2016a) and applied models both for igneous systems (e.g. Dufek and Bachmann, 2010) as well as engineering problems (e.g., Ho, 2005). The relevant physical scales in this limit can again be assessed using the generalised segregation and compaction numbers (58) as well as the segregation/compaction length (57). The small magnitude of δ s s/c shown in Fig. 6 indicates that crystal segregation speed is small compared to the convective flow of the carrier liquid. This insight is valid at least for systems much larger than the granular scale, 0 d 0 , as discussed before in the context of convecting magma chambers (Martin and Nokes, 1988;Brandeis and Jaupart, 1986;Rudman, 1992;Bergantz and Ni, 1999). Magma transport at high melt fractions hence is likely dominated by collective flow of partially crystalline magma with minor melt-crystal segregation.
The mechanical volume diffusion term in (59d) introduces a mechanical Peclet number based on the diffusive speed scale u s diff = K s φ,0 φ s 0 / 0 , with the mechanical diffusivity scaling as, K s φ,0 ∼ u s segr d 0 : The Peclet number thus scales with the cubed ratio of system to granular length scales, which means that phase diffusion is most consequential in systems not too far removed from the granular scale.
We have shown that the canonical limits of two-phase porous and suspension flow emerge as special cases of our multi-phase equations in the low-and high-melt fraction limits for incompressible solid-liquid systems with high phase viscosity contrast. Using the generalised segregation/compaction length scale, we have proposed two-phase regime boundaries delimiting the porous and suspension flow regimes as well as the intermediate mush regime. However, we have so far only touched on solid-liquid aggregates of melt and crystal/rock phases. Upper crustal magma processing, volcanic activity, and the genesis of mineral resources are, to a large degree, consequences of the exsolution of a third, volatile-rich, liquid to gas phase during decompression and cooling of magma in the shallow crust.

Three-phase flows in degassing magma
The solubility of volatiles, most notably H 2 O, CO 2 , and SO 2 , in igneous melts is a strong function of pressure but also depends on temperature and composition (Papale, 1999;Papale et al., 2006). As magma rises through the upper crust or cools in shallow magma bodies, volatiles exsolve from the melt to form different magmatic volatile phases: supercritical fluids, brines, hydrothermal liquids, and gaseous vapours. For simplicity, we collectively refer to that class of phases as volatile vapours (phase properties denoted as ( · ) v ) and will characterise them as generally very low-viscosity, low-density, compressible materials.
The challenge of building three-phase models in our framework lies in choosing closures that characterise the multi-phase interactions between vapour, liquid and solid phases. Prescribing appropriate transmission functions to represent local-scale phase topologies in igneous three-phase aggregate can be difficult. It is beyond the scope of this contribution to discuss even the leadingorder aspects of equilibrium textures between melt, crystals and vapour bubbles. Instead, we will demonstrate how to construct a three-phase model in our framework by giving an example. For the three-phase example, we use the same solid and liquid phase parameters as above and add in a third phase. We assume that third phase to be a vapour that preferentially wets the liquid rather than the solid. As a consequence, the vapour should have a higher percolation threshold than the liquid in a solid matrix, and the solid should have a higher disaggregation threshold in a solid-vapour than in a solid-liquid mixture. If sufficient liquid is present, the vapour phase should form bubbles up to very high vapour fractions where we assume a connected vapour phase will form. These assumptions form the basis for our choice of transmission functions. Figure 7(a)-(c) shows the intra-phase transmissions of our three-phase example plotted over the three-phase space. The transmission express, where a contiguous solid or a connected melt or vapour phase are present. In the solid transmission, the disaggregation threshold is shifted towards lower solid fractions along the solid-vapour axis compared to the solid-liquid axis. This comes as a consequence of assuming the vapour phase is less effectively wetting grain boundaries than the liquid (i.e. higher dihedral angle between vapour and solid than liquid and solid). The liquid transmission prescribes the collapse of liquid films between vapour bubbles by a drop-off above φ v ≈ 0.8. On the liquid-solid axis, we retain the low liquid percolation threshold at 2 % melt fraction from the two-phase example. For the vapour phase, we prescribe a percolation threshold in a solid matrix at 10 % vesicularity. The shift from disconnected bubbles to a connected vapour phase along the vapour-liquid axis is fixed at 0.8 vapour fraction, congruous with the choice of liquid disconnection at that point. The fitting parameters for this example are listed in Table B1, and the corresponding inter-phase transmissions are shown in Fig. B1.
Based on these transmissions, we now calculate the flux and transfer coefficients for our a three-phase example. We again use a granular length scale of 3 mm, and assume a viscosity of 10 −5 Pas and density of 100 kg/m 3 for the vapour phase. We show effective phase viscosities, phase segregation, and phase compaction coefficients in B2. We use the coefficients to quantify some pertinent metrics regarding the leading-order three-phase dynamics resulting from the chosen closures. The effective aggregate viscosity in Fig. 7(d) shows a region of high strength where the solid remains contiguous (teal); a second region of lower strength (beige), and a small third regime of lowest strength (brown) mark where the melt and vapour assume the role of carrier phase, respectively. In the latter regime, rock fragments and/or melt droplets are dynamically suspended in a gas phase as for example the case in products of explosive volcanic eruptions (ash clouds, pyroclastic density currents, etc.). Clearly, in this regime turbulent flow is to be expected. Our model is not intended for this regime and thus we will not further comment on it, nor would we recommend our framework be used in that limit.
The vapour segregation speed scale, u v segr (Fig. 7(e)) shows a region of very low mobility (< m/yr) in the solid limit, where vapour is present in mostly isolated vesicles and porous flow is limited. Much of the phase space at moderate to high melt fractions is characterised by an intermediate vapour segregation speed (m/yr-m/hr), where vapour bubbles segregate by displacing the melt and suspended crystals around them. At moderate to low melt fractions (< 0.2) and vapour fractions above the percolation threshold, the vapour phase is allowed to form a connected topology in a solid matrix causing the vapour segregation speed to pick up significantly to order m/s. vapour segregation in this region is no longer hindered by the viscosity of the carrier melt but rather by its own, much lower viscosity. Such a transition from vapour bubbles to an interconnected network as a function of both the vesicularity and crystallinity has been observed in analogue experiments (e.g., Colombier et al., 2018). The transition, if further corroborated, will clearly be a pivotal factor governing degassing of subvolcanic magma reservoirs and will need to be considered when interpreting observed changes in eruption styles, gas flux, or ground deformation at volcanic centers.
The length scale of vapour segregation through the compacting magma ( Fig. 7(f)) varies from the granular scale to more than a thousand kilometers. The maximum values are found where a connected vapour phase is present in a contiguous solid matrix. In such systems, efficient degassing may occur over long distances without any significant matrix deformation. However, outside that solid-dominated region, our example prescribes vapour segregation bubbles migrating through and deforming the variably crystalline mush. The low s/c-length in the liquid-dominated region of phase space suggests that segregation by bubble migration is slow compared to convective magma flow driven by vapour buoyancy. The degassing of mushy to melt-rich magma would therefore be facilitated by convective flow of a bubbly suspension rather than the relative migration of bubbles through the magma. Again, this interpretation relies on the assumption of bubbles much smaller than the system length scale.
Similarly to the two-phase limits above, we can us the velocity and pressure reference weights (see Fig. B3) along with the analysis of segregation/compaction lengths to identify different physical limits for which to introduce reduced equations. In Appendix B, we show such equations for the two limits of three-phase porous flow of connected melt and vapour phases through a contiguous rock matrix and for three-phase suspension of disaggregated crystals and bubbles in a carrier melt.
To sum up our considerations, we have drawn up a tentative regime diagram marking the regions of three-phase space where different styles of magma degassing may be expected (at least in our present example system). A first regime (1), sees the onset of degassing by porous flow in a solid-dominated matrix. As for two-phase systems, this regime is marked by a strong increase in segregation/compaction length, meaning that an addition of vapour increases both the segregation speed as well as the length scale of pressure-communicating porous medium. In the second regime (2) the vapour segregation speed through a disaggregating, melt-poor granular medium increases continuously. Whereas degassing is effective in this regime, the matrix becomes more deformable as the vapour fraction increases, until it may fragment, fluidise, and/or become entrained in the fast flowing vapour. The third regime (3) is one of bubbly flow in moderately crystalline to crystal-rich mushes. In this regime, degassing rates become rather sluggish and the compaction length drops gradually down to the granular scale as the mush becomes more mobile. The fourth regime (4) is characterised by crystal-poor bubbly suspension flow, where the s/c-length remains near the granular scale, and convective mobility exceeds segregation mobility of vapour bubbles for systems much larger than the mean bubble size. The final regime (5) is the aforementioned gas-carried suspension or ash cloud regime, where our model framework is no longer applicable.

Summary and Conclusions
In summary, we have derived a mixture theory framework for igneous reactive transport processes bridging the regimes of porous and suspension flows and extending existing two-phase models to n-phase generality. Our derivation draws heavily on both concepts from Rational Thermodynamics and procedures from Non-equilibrium Thermodynamics, which implies that we forgo a first-principles-based connection between local-scale phase interactions and system-scale processes. Instead, we base the model on spatially averaged conservation equations for interpenetrating phase continuum fields and assume that energies and entropies are meaningfully defined for the mixture of phases at the continuum scale.
Our model is formulated entirely at the macroscopic scale. An important consequence is that we have to choose constitutive relations that express system-scale mechanical and thermodynamic processes arising from local-scale phase interactions. We adhere to an internally consistent set of axiomatic principles and limiting assumptions to formulate constitutive relations for fluxes within and transfers between phases. The result is a set of linear, decoupled, isotropic relations that are admissible but non-unique. Our choices are guided by the motivation to find simple functional forms that satisfy the entropy inequality and align with basic empirical and conceptual thought on how multi-phase aggregates respond to two types of forces: spatial gradients and phase differences in thermodynamic variables.
The constitutive relations are chosen such that fluxes and transfers are proportional to their conjugate forces. The proportionality is set by material response coefficients, for which we propose a set of internally consistent but non-unique phenomenological closures, which require problemspecific calibration against field observations and laboratory data. The proposed closures are based on the broad assumption that fluxes and transfers are facilitated by the diffusion of heat, chemical species, momentum and partial volume within phases and between control volumes (fluxes), and between phases within a control volume (transfers). Applying this model to a specific natural or engineering context requires specification of pure-phase properties of the materials involved, as well as a set of transmission functions. The transmission functions capture both the iconnectivity of a phase within itself, as well as its connectivity with other phases. The latter aspect is particularly important in the context of vapour-bearing three-phase flows and allows distinguishing between bubbly and thin-film flows at the same relative vapour fraction. Together, the transmission functions and ratios of pure-phase properties are used to define flux and transfer permission functions that encapsulate the system-scale effects of local-scale material properties and phase topologies.
We have used the framework to construct examples of igneous two-and three-phase scenarios to discuss some ramifications for modelling the various stages of igneous processes from source to surface. With the proposed phenomenological closures, the general governing equations reduce to the classical limits of two-phase porous and suspension flow. We compare our closures to a variety of previous parameterizations including a Kozeny-Carman permeability and a hindered Stokes settling coefficient and to a variety of empirical and theoretical relations for effective viscosity (Krieger and Dougherty, 1959;Hirth and Kohlstedt, 2003;Costa et al., 2009;Takei and Holtzman, 2009;Rudge, 2018b) and mechanical diffusivity (Segre et al., 2001). Whereas we have focused on viscous materials here, our model framework lends itself to extension into the visco-elastic and brittleplastic domain. The additional rheological complexity will become necessary for modelling reactive transport across the ductile-brittle transition in crustal rocks. Finally, we have demonstrated how our model framework can be applied to three-phase flows involving a vapour phase in addition to a liquid and solid phase and draw up a tentative regime diagram for different expected modes of magma degassing.
Our model formulation intentionally separates the general governing equations from problemspecific material response coefficients. We argue that model transparency is increased by clearly distinguishing model components that are based on fundamental principles from phenomenological components that require calibration for a specific context. In practice, the steps of writing a physically and thermodynamically consistent set of equations and applying it to a specific application context by choosing an internally consistent set of closures are often conflated. The chief danger lies in the temptation to use a reduced set of governing equations designed for a particular limiting case (e.g., porous flow) outside the appropriate range of conditions (e.g., beyond solid disaggregation). We have therefore kept the derivation of the n-phase general governing equations separate from the model choices necessary to apply them to specific igneous problems. By allowing to transparently and flexibly formulate igneous models from a more general starting point, it is our hope that this framework will help the community find process-based common ground between the wide range of igneous phenomena from source to surface, as well as to establish better connections to other natural and engineering concerns of similar physical description.

Author Contributions
This study was a collaborative effort. TK conceived the study and took the lead in deriving the model and writing the manuscript. JS provided theoretical background of multi-phase continuum models and their connection to the local scale, and contributed to model formulation and writing.

A Assembling the Final Governing Equations
A.1 Phase segregation equation Starting from phase momentum conservation (6c), we substitute the constitutive relations for momentum transfer (20c), momentum diffusion flux (26c), and external momentum source (34b) and write the result in a form that emphasizes the role of phase velocity differences that govern phase segregation, Rearranging the first two terms in parentheses by substituting P * ∇φ i = ∇φ i P * − φ i ∇P * , we obtain, A.2 Phase compaction equation Next, we substitute the constitutive relations for volume transfer (20d), volume flux (26d), into phase mass conservation (6a) and write the result in a form that emphasizes the role of phase pressure differences that govern phase compaction, Rearranging the first two terms in parentheses by substituting v * · ∇φ i = ∇· φ i v * − φ i ∇· v * , we obtain,

A.3 Evolution of phase fractions
To write an equation for the evolution of phase fractions, φ i , we rearrange phase mass conservation (6a), and substitute the constitutive relation for the volume flux (26d), Alternatively expressing transport of phase fractions by substituting v i = v * + ∆v i * , we obtain, where the material derivative transported on v * is,

A.4 Evolution of component concentrations
We find the final equation for the evolution of component concentrations by substituting the constitutive relation for component fluxes (26b),

A.5 Evolution of phase temperatures
Energy conservation couples together all fluxes, transfers and sources. The procedure of substituting all constitutive relations and grouping terms to arrive at the final evolution equation becomes an exercise in bookkeeping. We begin with the initial Lagrangian form for conservation of total specific energy per phase (6d), into which we substitute the definition of total energy (8), In (9) internal energy is a function of specific entropy, volume, and chemical potentials. For the final form of the energy equation, it is useful to instead express it in terms of temperature, pressure and component concentrations (e.g., Rudge et al., 2011), with α i the phase thermal expansivity. We now substitute (A9) into (A8) to write a conservation equation for sensible heat coupled to the conservation of momentum, phase mass and component mass, Once more, we have sorted terms according to their categories of fluxes, transfers, and sources.
Next, we substitute conservation of momentum (6c), component mass (6b), and phase mass (6a) for the evolution terms on the right hand side of (A10): We have now assembled a temperature equation featuring all fluxes, transfers and sources, except for entropy terms, which for now are contained within energy terms. We now proceed to substitute the chosen constitutive relations, beginning with external sources (34), noting that potential energy terms cancel out. Next, we substitute constitutive relations for fluxes, (24) and (26), and transfers, (17) and (20), Note that two terms ±P i v i · ∇φ i arising from substituting fluxes and transfers have canceled out, and that once again a nonlinear disequilibrium term, ∆P i * ∆v i * · ∇φ i , was dropped. We further simplify (A14) by neglecting dissipation terms of component and volume fluxes, and component transfer, and introduce the latent heat of phase change, L i = T i s i , to obtain the final form, The final tally of processes contributing to temperature evolution are (from left to right), thermal diffusion and equilibration, dissipation of deviatoric, segregation, and compaction flow, latent heat of reaction, adiabatic heating, and radiogenic heating. Table B1 lists material and fitting parameters for the three-phase example discussed in the main text. All intra-and inter-phase transmissions resulting from these parameters are shown in Fig.  B1, and the corresponding effective phase viscosities, segregation and compaction coefficients in Fig. B2. Figure B3 finally shows the coefficient-based weights for reference velocity and pressure along with the pertinent segregation/compaction lengths. Based on the tentative regimes we identify from this example (see 8), and in analogy to the limiting cases of two-phase porous and suspension flow, we write governing equations for the two special limits of three-phase flow in solid-and liquid-dominated igneous aggregates. Although the complex equations of state for volatile vapor phases have received much attention (e.g., Sverjensky et al., 2014), we choose a linear pressure-dependent equation of state for simplicity: ρ v = ρ v 0 (1 + β v P v ), with β v [1/Pa] the isothermal vapor compressibility.

B.1 Three-phase Porous Flow
In the limit of three-phase porous flow (connected liquid and gas phases in a contiguous solid matrix), Fig. B3 shows that the solid velocity and vapour pressure assume the respective roles as reference velocitya and pressure. Hence, the mechanical governing equations can be reduced to, The addition of compressibility terms leads to a new dimensionless number that characterises the importance of compressibility relative to buoyancy-related pressure differences, R i compr = β i 0 ∆ρ ik g 0 δ ik s/c . For example, with a vapor compressibility of 10 −6 Pa −1 , and a vapor-solid density contrast of 2400 kg/m 3 , compressibility effects become significant if the s/c-length increases above 42 m.

B.2 Three-phase Suspension Flow
In the limit of three-phase suspension flow (disconnected solid particles and vapor bubbles suspended in a connected carrier melt), Fig. B3 shows that the liquid velocity assumes the reference role. Pressure differences are neglected for the same reasons as in the two-phase suspension limit. Hence, the mechanical governing equations can be reduced to,   Figure B1: Transmission functions for the example igneous three-phase system referenced in the main text. Brown colors denote high connectivity, while teal denotes disconnected phases. Figure B2: Momentum flux, segregation, and compaction coefficients for the example igneous threephase system referenced in the main text. Brown colors denote high phase mobility, while teal denote high resistance to deformation. Figure B3: Coefficient-based weights for reference velocity and pressure, and the liquid-in-solid, vaporin-solid, and vapor-in-liquid segregation-compaction lengths for the example igneous three-phase system referenced in the main text. Brown colors denote high phase mobility, while teal denote high resistance to deformation.