2.1. Fundamental assumptions and notation

We consider a continuum model of a mixture composed of two phases of a single chemical constituent. For our purposes, we suppose that the two phases are solid and liquid, labelled by superscripts $\varepsilon \in \{s, l\}$. Each of the two phases is a continuum thermodynamic system. The two thermodynamic continua are immiscible but interpenetrating (occupy the same space), and interact with each other in the form of mass, momentum and energy exchanges. We assume that all exchanges are pure, i.e. their sum vanishes. As such we neglect the effects of surface tension. Surface tension can be incorporated into the two-phase model (see e.g. Bercovici et al. Reference Bercovici, Ricard and Schubert2001) at the expense of additional complexity, which we feel is unwarranted at this stage of the model development and given that the interface properties of iron at Earth's core conditions (hundreds of GPa and several thousand K) are poorly known. We also assume that the two phases do not interact chemically. In other words, thermodynamic quantities associated with only one phase depend exclusively on the variables of this phase. For example, the Gibbs free energy of one phase is independent of the amount of the other phase present. This assumption also means that entropy of mixing is neglected.

At each point in the continuum the phases are characterised by their system-scale velocity $\boldsymbol {u}^\varepsilon$, pressure $p^\varepsilon$, temperature $T^\varepsilon$, etc. The relative amount of each phase at any point is measured by volume or mass fraction. These can be conceptually understood by considering a small finite control volume around a point $\boldsymbol {x}$ in space, of total mass $\mathcal {M}$ and volume $\mathcal {V}$. The total mass $\mathcal {M}$ is the sum of masses of each phase and, similarly, total volume $\mathcal {V}$ is the sum of volumes occupied by individual phases:

(2.1a,b)\begin{equation} \mathcal{M} = \mathcal{M}^s + \mathcal{M}^l, \quad \mathcal{V} = \mathcal{V}^s+\mathcal{V}^l.\end{equation}

Mass and volume fractions can thus be defined as

(2.2a,b)$$\begin{gather} \psi^\varepsilon = \frac{\mathcal{M}^\varepsilon}{\mathcal{M}}, \quad \phi^\varepsilon = \frac{\mathcal{V}^\varepsilon}{\mathcal{V}}; \end{gather}$$

(2.3a,b)$$\begin{gather}\psi^s + \psi^l = 1 , \quad \phi^s +\phi^l = 1. \end{gather}$$

Strictly speaking, at a microscopic level, any point $\boldsymbol {x}$ in space can only be occupied by either a solid phase or a liquid phase and thus volume/mass fractions should only take values 0 or 1. However, in the continuum approach, volume and mass fractions are interpreted as statistical averages, that is, $\phi ^s$ represents the fraction of a small volume surrounding $\boldsymbol {x}$ that is occupied by solid phase. Volume and mass fractions are thus continuous and differentiable functions of $\boldsymbol {x}$. This is the standard approach in continuum modelling of two-phase flows (e.g. Drew Reference Drew1983; Roberts & Loper Reference Roberts and Loper1987; Bercovici et al. Reference Bercovici, Ricard and Schubert2001).

The specific density $\rho ^\varepsilon$ is defined by the mass of phase $\varepsilon$ divided by the volume of phase $\varepsilon$, and is the inverse of the specific volume $V^\varepsilon$:

(2.4)\begin{equation} \rho^\varepsilon = \frac{\mathcal{M}^\varepsilon}{\mathcal{V}^\varepsilon} = \frac{1}{V^\varepsilon}.\end{equation}

In addition to quantities relating to individual phases, it is also useful to define quantities relating to the mixture. Throughout the paper we denote mixture quantities with an overbar. The specific density and specific volume of the mixture are

(2.5a,b)\begin{equation} \bar{\rho} = \frac{\mathcal{M}}{\mathcal{V}} = \phi^s \rho^s+\phi^l \rho^l, \quad \bar{V} = \frac{1}{\bar{\rho}} = \psi^s V^s + \psi^l V^l.\end{equation}

It is useful also to note the relation between mass fraction $\psi ^\varepsilon$ and volume fraction $\phi ^s$:

(2.6)\begin{equation} \psi^\varepsilon = \frac{\phi^\varepsilon \rho^\varepsilon}{\bar{\rho}}.\end{equation}

The velocity of the mixture $\boldsymbol {\bar {u}}$ is the centre-of-mass (barycentric) velocity defined by

(2.7)\begin{equation} \overline{\rho \boldsymbol{u}} = \phi^s \rho^s \boldsymbol{u}^s + \phi^l \rho^l \boldsymbol{u}^l.\end{equation}

Equivalently, $\boldsymbol {\bar {u}} = \overline {\rho \boldsymbol {u}}/\bar {\rho } = \psi ^s \boldsymbol {u}^s + \psi ^l \boldsymbol {u}^l$. Following Roberts & Loper (Reference Roberts and Loper1987) we suppose that the internal energy density of the mixture can be obtained by adding together those of the two phases separately:

(2.8)\begin{equation} \overline{\rho E} = \phi^s \rho^s E^s + \phi^l \rho^l E^l,\end{equation}

where $E^\varepsilon$ is the internal energy per unit mass of phase $\varepsilon$. This ignores interfacial surface tension. We assume that the entropy density of the mixture is similarly additive:

(2.9)\begin{equation} \overline{\rho S} = \phi^s \rho^s S^s + \phi^l \rho^l S^l,\end{equation}

where $S^\varepsilon$ is the specific entropy of phase $\varepsilon$. This form implicitly ignores entropy of mixing and interfacial surface tension.

We define material derivatives

(2.10a,b)\begin{equation} \frac{\mathrm{D}^\varepsilon }{\mathrm{D} t} \equiv \frac{\partial }{\partial t} + \boldsymbol{u}^\varepsilon \boldsymbol{\cdot} \boldsymbol{\nabla}, \quad \frac{\bar{\mathrm{D}} }{\bar{\mathrm{D}} t} \equiv \frac{\partial }{\partial t} + \boldsymbol{\bar{u}} \boldsymbol{\cdot} \boldsymbol{\nabla}, \end{equation}

where ${\mathrm {D}^\varepsilon }/{\mathrm {D} t}$ is the material derivative following the flow velocity of phase $\varepsilon$ and ${\bar {\mathrm {D}} }/{\bar {\mathrm {D}} t}$ is the material derivative following the barycentric velocity (2.7). We also introduce notation for phase difference in any quantity $A$:

(2.11)\begin{equation} \Delta A = A^s - A^l. \end{equation}

We derive all equations to maintain their Galilean and material invariance. Galilean invariance requires that the laws of motion are independent of the inertial frame of reference. Material invariance refers to the symmetry between the equations governing the two phases, whereby an interchange of labels results in the same equations.

2.2. Conservation of mass

The mass continuity equation for phase $\varepsilon$ is

(2.12)\begin{equation} \frac{\partial (\phi^\varepsilon \rho^\varepsilon)}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^\varepsilon \rho^\varepsilon \boldsymbol{u}^\varepsilon ) = \varGamma^\varepsilon,\end{equation}

where $\varGamma ^\varepsilon$ is the rate of production of phase $\varepsilon$. The rate of production of one phase is equal and opposite to the rate of destruction of the other, i.e. $\varGamma ^s + \varGamma ^l = 0$. The total mass conservation law is obtained by summing continuity equations of individual phases together:

(2.13)\begin{equation} \frac{\partial \bar{\rho}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \overline{\rho \boldsymbol{u}} ) = 0.\end{equation}

Using (2.6) and (2.13) in (2.12) we obtain an evolution equation for solid mass fraction $\psi ^s$:

(2.14)\begin{equation} \bar{\rho} \frac{\bar{\mathrm{D}} \psi^s}{\bar{\mathrm{D}} t} + \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{j} = \varGamma^s, \quad \boldsymbol{j} = \bar{\rho} \psi^s (1-\psi^s)\Delta \boldsymbol{u} = \frac{\phi^s \rho^s \phi^l \rho^l}{\bar{\rho}} \Delta \boldsymbol{u}, \end{equation}

where $\boldsymbol {j}$ is the solid flux (Loper & Roberts Reference Loper and Roberts1978).

Note also that the phase mass continuity equation (2.12) can be written as

(2.15)\begin{equation} \frac{\mathrm{D}^\varepsilon \phi^\varepsilon }{\mathrm{D} t} =\frac{\varGamma^\varepsilon}{\rho^\varepsilon} -\phi^\varepsilon \zeta^\varepsilon , \quad \text{where } \zeta^\varepsilon = \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}^\varepsilon + \frac{1}{\rho^\varepsilon} \frac{\mathrm{D}^\varepsilon \rho^\varepsilon}{\mathrm{D} t} .\end{equation}

Equation (2.15) states that the volume fraction $\phi ^\varepsilon$ can increase as a result of production of the phase ($\varGamma ^\varepsilon > 0$), convergence of the phase due to flow ($\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {u}^\varepsilon <0$) or a decrease in specific density (${\mathrm {D}^\varepsilon \rho ^\varepsilon }/{\mathrm {D} t} < 0$; equivalently, an increase in specific volume ${\mathrm {D}^\varepsilon V^\varepsilon }/{\mathrm {D} t} > 0$). Solid compaction is the process by which a solid matrix progressively loses its porosity due to the effects of pressure from loading – i.e. compaction decreases $\phi ^l$ and increases $\phi ^s$. We can generalise this terminology to say that compaction of phase $\varepsilon$ is a process which increases $\phi ^\varepsilon$. Since compaction is a process that is treated as separate from melting or freezing, the quantity $\zeta ^\varepsilon$ thus denotes compaction ($\zeta ^\varepsilon < 0$) and dilation ($\zeta ^\varepsilon > 0$) of phase $\varepsilon$.

2.3. Conservation of momentum

The general form of the equation for the rate of change of momentum is

(2.16)\begin{equation} \frac{\partial }{\partial t}( \phi^\varepsilon \rho^\varepsilon \boldsymbol{u}^\varepsilon ) + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^\varepsilon \rho^\varepsilon \boldsymbol{u}^\varepsilon \boldsymbol{u}^\varepsilon ) = \boldsymbol{F}^\varepsilon ,\end{equation}

where $\boldsymbol {F}^\varepsilon$ stands for the totality of forces through which momentum is gained/redistributed within phase $\varepsilon$. The momentum transfer term consists of contributions due to internal stress (due to pressure and viscous forces within the phase), external forces (e.g. gravity) and momentum exchange due to interaction between phases:

(2.17)\begin{equation} \boldsymbol{F}^\varepsilon = \boldsymbol{F}^\varepsilon_{{internal}} + \boldsymbol{F}^\varepsilon_{{external}} + \boldsymbol{F}^\varepsilon_{{interaction}} .\end{equation}

The internal force arises due to divergence of the total stress tensor $\boldsymbol {\varPi }^\varepsilon$:

(2.18)\begin{equation} \boldsymbol{F}_{internal}^\varepsilon = \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^\varepsilon \boldsymbol{\varPi}^\varepsilon ) ,\end{equation}

where it has been assumed that the stress tensors of the individual phases are additive. The total stress tensor for phase $\varepsilon$ can be split into isotropic and deviatoric components:

(2.19)\begin{equation} \boldsymbol{\varPi}^\varepsilon =- p^\varepsilon \boldsymbol{ I } + \boldsymbol{\sigma}^\varepsilon,\end{equation}

where $p^\varepsilon$ is the mechanical pressure (minus one-third trace of the stress tensor, $p^\varepsilon = - \frac {1}{3} \mathrm {Tr}( \boldsymbol {\varPi }^\varepsilon )$), $\boldsymbol { I }$ is the identity matrix and $\boldsymbol {\sigma }^\varepsilon$ is the deviatoric stress tensor. Force densities on each phase due to the external gravitational field $\boldsymbol { g }$ are

(2.20)\begin{equation} \boldsymbol{F}^\varepsilon_{{external}} = \phi^\varepsilon \rho^\varepsilon \boldsymbol{ g }.\end{equation}

The interaction force $\boldsymbol {F}^\varepsilon _{{interaction}}$ results from forces acting at the interface between the phases. Given the complexity of the interface between the phases, and lack of knowledge of its precise location and orientation in the continuum approach, the interfacial forces are difficult to quantify and have been subject of much debate in the two-phase literature. The interaction force is commonly split into three components: momentum exchange due to transfer of mass between the phases during melting and solidification; a ‘buoyant’ force each phase experiences due to the pressure exerted on its boundary by the other phases which surround it; and a frictional force that phases experience as they move past one another (Drew Reference Drew1983). The momentum exchange due to melting and freezing depends on the velocity at which the solid–liquid interface propagates during the phase change. Similarly, the buoyant force between the phases depends on the pressure distribution on solid–liquid interfaces. Neither the ‘interface velocity’ nor the ‘interface pressure’ are known (or even knowable) quantities in the continuum framework, where even the concept of the interface is blurred (in the continuum framework interfaces are not discrete entities; the interface as well as the two phases are continuous quantities that exist at all points in the domain), and have to be parametrised in terms of phase velocities $\boldsymbol {u}^\varepsilon$ and phase pressures $p^\varepsilon$. From these considerations, the interaction force can be expressed as (see e.g. Ishii & Hibiki (Reference Ishii and Hibiki2010), chapter 9)

(2.21)\begin{equation} \boldsymbol{F}^\varepsilon_{{interaction}} = \varGamma^\varepsilon \boldsymbol{\tilde{u}}^\varepsilon + \tilde{p}^\varepsilon \boldsymbol{\nabla} \phi^\varepsilon + \boldsymbol{D}^\varepsilon , \end{equation}

where $\boldsymbol {\tilde {u}}^\varepsilon$ is the interface velocity on the side of phase $\varepsilon, \tilde {p}^\varepsilon$ is the interface pressure on the side of phase $\varepsilon\ {\rm and}\ \boldsymbol {D}^\varepsilon$ is the generalised drag force that is equal and opposite in each phase ($\boldsymbol {D}^s =-\boldsymbol {D}^l$). In the absence of surface tension the interaction forces are equal and opposite, $\boldsymbol {F}^s_{{interaction}} = -\boldsymbol {F}^l_{{interaction}}$, and interface velocity and pressure are continuous across the phase interfaces:

(2.22a,b)\begin{equation} \boldsymbol{\tilde{u}}^s = \boldsymbol{\tilde{u}}^l = \boldsymbol{\tilde{u}}, \quad \tilde{p}^s = \tilde{p}^l = \tilde{p}.\end{equation}

Utilising the phase mass conservation equation (2.12), and relations (2.17)–(2.22a,b), the momentum equation (2.16) becomes

(2.23) \begin{align} \phi^\varepsilon \rho^\varepsilon \left(\frac{\partial \boldsymbol{u}^\varepsilon}{\partial t} + \boldsymbol{u}^\varepsilon \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}^\varepsilon \right) &=- \phi^\varepsilon \boldsymbol{\nabla} p^\varepsilon + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^\varepsilon \boldsymbol{\sigma}^\varepsilon ) + \phi^\varepsilon \rho^\varepsilon \boldsymbol{ g } + \varGamma^\varepsilon ( \boldsymbol{\tilde{u}} - \boldsymbol{u}^\varepsilon )\notag\\ &\quad + (\tilde{p} - p^\varepsilon ) \boldsymbol{\nabla} \phi^\varepsilon + \boldsymbol{D}^\varepsilon. \end{align}

Following Roberts & Loper (Reference Roberts and Loper1987) and Bercovici et al. (Reference Bercovici, Ricard and Schubert2001) we postulate the following constitutive relations for $\boldsymbol {\tilde {u}}$ and $\tilde {p}$:

(2.24a–c)$$\begin{gather} \boldsymbol{u}^s - \boldsymbol{\tilde{u}} = \gamma^s \Delta \boldsymbol{u}, \quad \boldsymbol{u}^l - \boldsymbol{\tilde{u}} =-\gamma^l \Delta \boldsymbol{u}, \quad \gamma^s +\gamma^l = 1; \end{gather}$$

(2.25a–c)$$\begin{gather}p^s - \tilde{p} = \omega^s \Delta p , \quad p^l - \tilde{p} =-\omega^l \Delta p, \quad \omega^s + \omega^l = 1. \end{gather}$$

Coefficients $\omega ^\varepsilon$ control the partitioning of the interface surface force between the two phases and represent the fraction of the volume-averaged surface force exerted on each phase (Bercovici & Ricard Reference Bercovici and Ricard2003). Similarly, coefficients $\gamma ^\varepsilon$ control the partitioning of the interface velocity $\boldsymbol {\tilde {u}}$. On application of (2.24a,b) and (2.25a,b) in (2.23), momentum equations for solid and liquid phases become

(2.26a) $$\begin{gather} \phi^s \rho^s \left( \frac{\partial \boldsymbol{u}^s}{\partial t} {\,+\,} \boldsymbol{u}^s \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}^s \right) {\,=\,}- \phi^s \boldsymbol{\nabla} p^s {\,+\,} \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^s \boldsymbol{\sigma}^s ) {\,+\,} \phi^s \rho^s \boldsymbol{ g } {\,-\,} \gamma^s \varGamma^s \Delta \boldsymbol{u} {\,-\,} \omega^s \Delta p \boldsymbol{\nabla} \phi^s {\,+\,} \boldsymbol{D}^s, \end{gather}$$

(2.26b)$$\begin{gather}\phi^l \rho^l \left( \frac{\partial \boldsymbol{u}^l}{\partial t} + \boldsymbol{u}^l \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}^l \right) =- \phi^l \boldsymbol{\nabla} p^l + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^l \boldsymbol{\sigma}^l ) + \phi^l \rho^l \boldsymbol{ g } + \gamma^l \varGamma^l \Delta \boldsymbol{u} + \omega^l \Delta p \boldsymbol{\nabla} \phi^l + \boldsymbol{D}^l. \end{gather}$$

Taking the dot product of (2.26a) with $\boldsymbol {u}^s$ and the dot product of (2.26b) with $\boldsymbol {u}^l$, we obtain the equations for the evolution of kinetic energy of the two phases:

(2.27a)\begin{align} &\frac{\partial \mathcal{K}^s}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \mathcal{K}^s \boldsymbol{u}^s ) =- \phi^s \boldsymbol{u}^s \boldsymbol{\cdot} \boldsymbol{\nabla} p^s + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^s \boldsymbol{u}^s \boldsymbol{\cdot} \boldsymbol{\sigma}^s ) - \phi^s \boldsymbol{\sigma}^s : \boldsymbol{\nabla} \boldsymbol{u}^s\nonumber\\ &\quad + \phi^s \rho^s \boldsymbol{u}^s \boldsymbol{\cdot} \boldsymbol{ g } +\varGamma^s \boldsymbol{u}^s \boldsymbol{\cdot} \left( \frac{1}{2} \boldsymbol{u}^s - \gamma^s \Delta \boldsymbol{u} \right) - \Delta p \omega^s \boldsymbol{u}^s \boldsymbol{\cdot} \boldsymbol{\nabla} \phi^s + \boldsymbol{u}^s \boldsymbol{\cdot} \boldsymbol{D}^s, \end{align}

(2.27b)\begin{align} &\frac{\partial \mathcal{K}^l}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \mathcal{K}^l \boldsymbol{u}^l ) =- \phi^l \boldsymbol{u}^l \boldsymbol{\cdot} \boldsymbol{\nabla} p^l + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^l \boldsymbol{u}^l \boldsymbol{\cdot} \boldsymbol{\sigma}^l ) - \phi^l \boldsymbol{\sigma}^l : \boldsymbol{\nabla} \boldsymbol{u}^l\nonumber\\ &\quad + \phi^l \rho^l \boldsymbol{u}^l \boldsymbol{\cdot} \boldsymbol{ g } + \varGamma^l \boldsymbol{u}^l \boldsymbol{\cdot} \left( \frac{1}{2} \boldsymbol{u}^l + \gamma^l\Delta \boldsymbol{u} \right) + \Delta p \omega^l \boldsymbol{u}^l \boldsymbol{\cdot} \boldsymbol{\nabla} \phi^l + \boldsymbol{u}^l \boldsymbol{\cdot}\boldsymbol{D}^l, \end{align}

where $\mathcal {K}^\varepsilon = \frac {1}{2} \phi ^\varepsilon \rho ^\varepsilon | \boldsymbol {u}^\varepsilon |^2$ is the kinetic energy density of phase $\varepsilon$ and $\boldsymbol {\sigma }^\varepsilon : \boldsymbol {\nabla } \boldsymbol {u}^\varepsilon = \sigma ^\varepsilon _{ij} \partial {u_i^\varepsilon }/\partial {x_j}$.

In general, the gravity vector can be expressed in terms of the gradient of gravitational potential $\boldsymbol { g } = -\boldsymbol {\nabla } \varPhi$. Thus, with use of the continuity equation (2.12), and assuming that the gravitational potential is time-independent $\partial _{t} \varPhi = 0$, we can express the gravitational term in the kinetic energy equation as follows:

(2.28)\begin{align} \phi^\varepsilon \rho^\varepsilon \boldsymbol{u}^\varepsilon \boldsymbol{\cdot} \boldsymbol{ g } &=-\phi^\varepsilon \rho^\varepsilon \boldsymbol{u}^\varepsilon \boldsymbol{\cdot} \boldsymbol{\nabla} \varPhi =-\boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^\varepsilon \rho^\varepsilon \boldsymbol{u}^\varepsilon \varPhi ) + \varPhi \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^\varepsilon \rho^\varepsilon \boldsymbol{u}^\varepsilon ) \nonumber\\ &=-\left( \frac{\partial \mathcal{P}^\varepsilon }{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \mathcal{P}^\varepsilon \boldsymbol{u}^\varepsilon ) \right) +\varGamma^\varepsilon \varPhi , \end{align}

where $\mathcal {P}^\varepsilon = \phi ^\varepsilon \rho ^\varepsilon \varPhi$ is the potential energy density of phase $\varepsilon$. Adding together the kinetic energy equations of the two phases (2.27a) and (2.27b) we obtain the equation governing the evolution of total mechanical (kinetic and potential) energy in the mixture:

(2.29)\begin{align} &\frac{\partial }{\partial t} ( \mathcal{K}^s+ \mathcal{P}^s + \mathcal{K}^l + \mathcal{P}^l ) + \boldsymbol{\nabla}\boldsymbol{\cdot} ( ( \mathcal{K}^s+ \mathcal{P}^s )\boldsymbol{u}^s+ ( \mathcal{K}^l + \mathcal{P}^l )\boldsymbol{u}^l - \phi^s \boldsymbol{u}^s \boldsymbol{\cdot} \boldsymbol{\sigma}^s - \phi^l \boldsymbol{u}^l \boldsymbol{\cdot} \boldsymbol{\sigma}^l ) \nonumber\\ &\quad =- \phi^s \boldsymbol{u}^s \boldsymbol{\cdot} \boldsymbol{\nabla} p^s - \phi^l \boldsymbol{u}^l \boldsymbol{\cdot} \boldsymbol{\nabla} p^l - \phi^s \boldsymbol{\sigma}^s \! : \! \boldsymbol{\nabla} \boldsymbol{u}^s - \phi^l \boldsymbol{\sigma}^l \! : \! \boldsymbol{\nabla} \boldsymbol{u}^l + \boldsymbol{D}^s \boldsymbol{\cdot} \Delta \boldsymbol{u}\nonumber\\ &\qquad - \frac{1}{2} ( \varGamma^s \gamma^s + \varGamma^l \gamma^l ) | \Delta \boldsymbol{u} |^2 - \Delta p \boldsymbol{u}^\omega \boldsymbol{\cdot} \boldsymbol{\nabla} \phi^s, \end{align}

where $\boldsymbol {u}^\omega = \omega ^s \boldsymbol {u}^s + \omega ^l \boldsymbol {u}^l$. Terms on the left-hand side of (2.29) represent changes in mechanical energy due to accumulation, advective flux and viscous deformation stresses; terms on the right-hand side represent work against the pressure gradient, viscous heating, frictional heating, contributions due to phase changes and an additional contribution to the usual pressure–work term due to the existence of pressure disequilibrium ($\Delta p \neq 0$).

2.4. Conservation of entropy and energy

The entropy equation is

(2.30)\begin{equation} \frac{\partial }{\partial t}(\phi^\varepsilon \rho^\varepsilon S^\varepsilon) + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^\varepsilon \rho^\varepsilon S^\varepsilon \boldsymbol{u}^\varepsilon + \boldsymbol{ k }^\varepsilon ) = \varSigma^\varepsilon,\end{equation}

where $\boldsymbol { k }^\varepsilon$ represents the entropy flux due to internal processes (e.g. diffusive heat flux) and $\varSigma ^\varepsilon$ represents sources of entropy due to dissipative processes (e.g. viscous heating) and interaction with other phases. The second law of thermodynamics requires that the sum of the entropies of all bodies taking part in the process is increased, which is expressed as

(2.31)\begin{equation} \frac{\mathrm{d} }{\mathrm{d} t} \int _V \overline{\rho S} \, \mathrm{d} V = \int _V \varSigma^s + \varSigma^l \, \mathrm{d} V \geq 0. \end{equation}

This serves as a powerful constraint in developing constitutive laws.

To obtain the evolution equation for internal energy of each phase we begin with the statement of the first law of thermodynamics applied to the individual phases. According to the first law, the change in the internal energy is equal to net energy added in the form of heat $\delta q^\varepsilon = T^\varepsilon \mathrm {d} S^\varepsilon$ minus energy spent in the form of thermodynamic work $\delta w^\varepsilon = -p^\varepsilon \mathrm {d} V^\varepsilon$:

(2.32)\begin{equation} \mathrm{d} E^\varepsilon = T^\varepsilon \mathrm{d} S^\varepsilon - p^\varepsilon \mathrm{d} (1/\rho^\varepsilon) .\end{equation}

We actually want a relation for the internal energy density $\phi ^\varepsilon \rho ^\varepsilon E^\varepsilon$. We obtain this by multiplying (2.32) by $\phi ^\varepsilon \rho ^\varepsilon$ and rearranging:

(2.33)\begin{equation} \mathrm{d} (\phi^\varepsilon \rho^\varepsilon E^\varepsilon) = T^\varepsilon \mathrm{d} (\phi^\varepsilon \rho^\varepsilon S^\varepsilon) - p^\varepsilon \mathrm{d} \phi^\varepsilon + G^\varepsilon \mathrm{d} (\phi^\varepsilon \rho^\varepsilon),\end{equation}

where

(2.34)\begin{equation} G^\varepsilon = E^\varepsilon - T^\varepsilon S^\varepsilon + \frac{p^\varepsilon }{\rho^\varepsilon}\end{equation}

is the specific Gibbs free energy of phase $\varepsilon$. (Note that the specific Gibbs free energy introduced in (2.34) is often denoted by $\mu$ (e.g. Roberts & Loper (Reference Roberts and Loper1987), their (2.16)). Here we use $G$ to denote specific Gibbs free energy and reserve $\mu ^\varepsilon$ for dynamic viscosity.) Using (2.33) we obtain a material derivative of energy density $\rho ^\varepsilon E^\varepsilon$:

(2.35)\begin{equation} \frac{\mathrm{D}^\varepsilon (\phi^\varepsilon \rho^\varepsilon E^\varepsilon)}{\mathrm{D} t} = T^\varepsilon \frac{\mathrm{D}^\varepsilon (\phi^\varepsilon \rho^\varepsilon S^\varepsilon)}{\mathrm{D} t} - p^\varepsilon \frac{\mathrm{D}^\varepsilon \phi^\varepsilon}{\mathrm{D} t} + G^\varepsilon \frac{\mathrm{D}^\varepsilon (\phi^\varepsilon \rho^\varepsilon)}{\mathrm{D} t}.\end{equation}

Using (2.12), (2.30) and (2.34) in (2.35) yields (after a few vector identities) the internal energy equation:

(2.36) \begin{align} \frac{\partial (\phi^\varepsilon \rho^\varepsilon E^\varepsilon)}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} [ ( \phi^\varepsilon \rho^\varepsilon E^\varepsilon + \phi^\varepsilon p^\varepsilon )\boldsymbol{u}^\varepsilon ] &=- \boldsymbol{\nabla}\boldsymbol{\cdot} (T^\varepsilon \boldsymbol{ k }^\varepsilon) + \boldsymbol{ k }^\varepsilon \boldsymbol{\cdot} \boldsymbol{\nabla} T^\varepsilon + T^\varepsilon \varSigma^\varepsilon \nonumber\\ &\quad\, +\phi^\varepsilon \boldsymbol{u}^\varepsilon \boldsymbol{\cdot} \boldsymbol{\nabla} p^\varepsilon + G^\varepsilon \varGamma^\varepsilon - p^\varepsilon \frac{\partial \phi^\varepsilon}{\partial t} . \end{align}

The equation governing the total energy of the system is obtained by summing the solid and liquid internal energy equations (2.36) with the total kinetic energy equation (2.29),

(2.37)\begin{equation} \frac{\partial (\phi^s \rho^s \mathcal{E}^s + \phi^l \rho^l \mathcal{E}^l)}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \mathcal{F}^s + \mathcal{F}^l ) = \mathcal{J},\end{equation}

where $\mathcal {E}^\varepsilon$ is the total energy of phase $\varepsilon$ (per unit mass of phase $\varepsilon$), $\mathcal {E}^\varepsilon = ( E^\varepsilon + \frac {1}{2} | \boldsymbol {u}^\varepsilon |^2 + \varPhi ), \mathcal {F}^\varepsilon$ is the energy flux vector, ${\mathcal {F}^\varepsilon = ( \phi ^\varepsilon \rho ^\varepsilon \mathcal {E}^\varepsilon + \phi ^\varepsilon p^\varepsilon ) \boldsymbol {u}^\varepsilon - \phi ^\varepsilon \boldsymbol {u}^\varepsilon \boldsymbol {\cdot } \boldsymbol {\sigma }^\varepsilon + T^\varepsilon \boldsymbol { k }^\varepsilon }$ and $\mathcal {J}$ contains the sum of energy exchange terms arising due to the interaction between the phases:

(2.38)\begin{align} \mathcal{J} &=- \phi^s \boldsymbol{\sigma}^s : \boldsymbol{\nabla} \boldsymbol{u}^s - \phi^l \boldsymbol{\sigma}^l : \boldsymbol{\nabla} \boldsymbol{u}^l + \boldsymbol{ k }^s \boldsymbol{\cdot} \boldsymbol{\nabla} T^s + \boldsymbol{ k }^l \boldsymbol{\cdot} \boldsymbol{\nabla} T^l + T^s \varSigma^s + T^l \varSigma^l \nonumber\\ &\quad + \boldsymbol{D}^s \boldsymbol{\cdot} \Delta \boldsymbol{u} + \Delta G \varGamma^s - \frac{1}{2} ( \varGamma^s\gamma^s + \varGamma^l\gamma^l ) | \Delta \boldsymbol{u} |^2 - \Delta p \left( \frac{\partial \phi^s}{\partial t} + \boldsymbol{u}^{\omega} \boldsymbol{\cdot} \boldsymbol{\nabla} \phi^s \right). \end{align}

The material derivative appearing in the final term of (2.38) can be obtained by multiplying equation (2.15) corresponding to the solid phase by $\omega ^s$, and that corresponding to the liquid phase by $-\omega ^l$, and adding the resulting two equations to yield

(2.39)\begin{equation} \left( \frac{\partial \phi^s}{\partial t} + \boldsymbol{u}^{\omega} \boldsymbol{\cdot} \boldsymbol{\nabla} \phi^s \right) = \varGamma^s \frac{\tilde{\rho}}{\rho^s \rho^l} -\omega^s \phi^s \zeta^s + \omega^l \phi^l \zeta^l, \quad \text{where}\ \tilde{\rho} = \omega^l\rho^s + \omega^s\rho^l.\end{equation}

In a closed system, conservation of total energy in the system requires that

(2.40)\begin{equation} \frac{\mathrm{d} }{\mathrm{d} t} \int _V \mathcal{E}^s + \mathcal{E}^l \, \mathrm{d} V = \int _V \mathcal{J} \, \mathrm{d} V = 0.\end{equation}

Interaction between the phases cannot change the overall energy of the mixture (only transfer the energy from one phase to another). This is in line with our (more stringent) assumption that the energy exchanges are pure (their sum vanishes), which requires that

(2.41)\begin{equation} \mathcal{J} = 0.\end{equation}

2.5. Thermodynamic relations

Before establishing thermodynamic relations for the two-phase system, we must address the subtle matter of pressure and in particular the difference between ‘thermodynamic’ and ‘mechanical’ definitions of pressure. Even for a compressible single-phase viscous fluid there is a difference between these two definitions of pressure. The mechanical pressure is defined by minus one third the trace of the stress tensor. The thermodynamic pressure appears in definitions of thermodynamic potentials and accounts for changes in internal energy with respect to volume in a compressible fluid. The difference between the two definitions depends on the compaction viscosity and the divergence of the velocity field, and disappears in static equilibrium (e.g. Landau & Lifshitz Reference Landau and Lifshitz1987).

So far we have introduced two mechanical pressures ($p^s, p^l$), one for each of the two phases, and a shared interface pressure ($\tilde {p}$). Following Ricard (Reference Ricard2007) (see also Rudge, Bercovici & Spiegelman Reference Rudge, Bercovici and Spiegelman2011), in what follows we identify the interface pressure as the appropriate thermodynamic pressure, which we redefine as $P=\tilde {p}$ . Thus, from (2.25a,b),

(2.42a,b)\begin{equation} p^s = P + \omega^s \Delta p , \quad p^l = P - \omega^l \Delta p .\end{equation}

The pressure of each phase is assumed to be composed of the sum of thermodynamic pressure $P$ and compaction pressure $\Delta p$, which is entirely mechanical in nature as it arises due to flow. All thermodynamic potentials for both phases are defined using the thermodynamic pressure $P$. All of the thermodynamic differential relations involving pressure are to be taken as differentials of thermodynamic pressure $P$, i.e.

(2.43)\begin{equation} \mathrm{d} p^s = \mathrm{d} p^l = \mathrm{d} P. \end{equation}

The derivation of thermodynamic relations for Gibbs free energy $G^\varepsilon$, specific density $\rho ^\varepsilon$ and specific entropy $S^\varepsilon$ is standard and has been relegated to Appendix B. These differential relations are as follows:

(2.44)$$\begin{gather} \mathrm{d} G^\varepsilon = \frac{1}{\rho^\varepsilon} \mathrm{d} P - S^\varepsilon \mathrm{d} T^\varepsilon, \end{gather}$$

(2.45)$$\begin{gather}\mathrm{d} \rho^\varepsilon = \rho^\varepsilon ( \beta^\varepsilon \mathrm{d} P -\alpha^\varepsilon \mathrm{d} T^\varepsilon ), \end{gather}$$

(2.46)$$\begin{gather}\mathrm{d} S^\varepsilon =-\frac{\alpha^\varepsilon}{\rho^\varepsilon} \mathrm{d} P +\frac{c_p^\varepsilon}{T^\varepsilon}\mathrm{d} T^\varepsilon , \end{gather}$$

where $\alpha ^\varepsilon$ is the thermal expansion coefficient, $\beta ^\varepsilon$ is the isothermal compression coefficient and $c_p^\varepsilon$ is the specific heat capacity at constant pressure of phase $\varepsilon$.

In light of the pressure decomposition (2.42), it is useful also to introduce the effective Gibbs free energy $G_*^\varepsilon$ as the Gibbs free energy defined using the thermodynamic pressure (Ricard Reference Ricard2007) (compare with $G^\varepsilon$ defined in (2.34)):

(2.47a–c)\begin{equation} G_*^\varepsilon = E^\varepsilon - T^\varepsilon S^\varepsilon + \frac{P }{\rho^\varepsilon}, \quad G^s = G_*^s + \frac{\omega^s \Delta p }{\rho^s} , \quad G^l = G_*^l - \frac{\omega^l \Delta p }{\rho^l}.\end{equation}

Thus the Gibbs free energy of each phase can be thought of as a sum of thermodynamic Gibbs free energy ($G_*^\varepsilon$) and a contribution due to compaction pressure $\Delta p$, which is entirely mechanical in nature. Between (2.47a–c) and (2.44) it immediately follows that

(2.48a)$$\begin{gather} \mathrm{d} G_*^s = \frac{1}{\rho^s} \mathrm{d} P - S^s \mathrm{d} T^s - \omega^s \Delta p \, \mathrm{d} \left( \frac{1}{\rho^s} \right), \end{gather}$$

(2.48b)$$\begin{gather}\mathrm{d} G_*^l = \frac{1}{\rho^l} \mathrm{d} P - S^l \mathrm{d} T^l + \omega^l \Delta p \, \mathrm{d} \left( \frac{1}{\rho^l} \right). \end{gather}$$

In the case of mechanical disequilibrium ($\Delta p \neq 0$), the chemical equilibrium is reached when the effective Gibbs free energies of the two phases are equal.

2.6. Constitutive relations

In this section we exploit the principles of energy conservation and entropy production (2.31) to derive constitutive relations for the closure terms in the governing equations. We need to provide expressions for the deviatoric stress tensor $\boldsymbol {\sigma }^\varepsilon$, heat flux vector $\boldsymbol { k }^\varepsilon$, frictional momentum exchange $\boldsymbol {D}^\varepsilon$, mass transfer term $\varGamma ^\varepsilon$, pressure difference between phases $\Delta p$, entropy production $\varSigma ^\varepsilon$ and coefficients $\gamma ^\varepsilon, \omega ^\varepsilon$.

We begin with a parametrisation of $\gamma ^\varepsilon$ since it is the most elusive one, and in the absence of physical estimates it has to be fixed based on theoretical reasoning. We recall that $\gamma ^\varepsilon$ is associated with the process of kinetic energy exchange between the phases during melting and solidification. We take a simple view that the kinetic energy exchange due to phase change is pure (kinetic energy of one phase is transferred wholly into the kinetic energy of the other phase, with no energy being converted into other forms of energy), and thus does not change the overall mechanical energy of the system. The terms involving $\gamma ^\varepsilon$ in the mechanical energy equation (2.29) should therefore vanish, which requires that

(2.49)\begin{equation} \gamma^s \varGamma^s + \gamma^l \varGamma^l = 0. \end{equation}

The constraint $\gamma ^s + \gamma ^l = 1$ (equation (2.24c)) fixes

(2.50)\begin{equation} \gamma^\varepsilon = \tfrac{1}{2}.\end{equation}

This parametrisation eliminates the $\gamma ^\varepsilon$-dependent energy source term in $\mathcal {J}$ (2.38). Alternative parametrisations are possible, for example setting $\gamma ^\varepsilon = \mathcal {H} ( \varGamma ^\varepsilon )$, where $\mathcal {H}$ is the Heaviside function (Roberts & Loper Reference Roberts and Loper1987), at the expense of increased mathematical complexity.

Conservation of total energy (2.40) requires that $\mathcal {J}$ vanishes. Thus on using (2.39) for ${\mathrm {D}^\omega \phi ^s}/{\mathrm {D} t}$ and (2.47a–c) for $\Delta G = \Delta G_* + \Delta p {\tilde {\rho }}/(\rho ^s \rho ^l)$ in (2.38) yields

(2.51)\begin{align} T^s \varSigma^s + T^l \varSigma^l &= \phi^s \boldsymbol{\sigma}^s : \boldsymbol{\nabla} \boldsymbol{u}^s + \phi^l \boldsymbol{\sigma}^l : \boldsymbol{\nabla} \boldsymbol{u}^l - \boldsymbol{ k }^s \boldsymbol{\cdot} \boldsymbol{\nabla} T^s - \boldsymbol{ k }^l \boldsymbol{\cdot} \boldsymbol{\nabla} T^l \nonumber\\ &\quad - \boldsymbol{D}^s \boldsymbol{\cdot} \Delta \boldsymbol{u} - \varGamma^s \Delta G_* + \Delta p ( -\omega^s \phi^s \zeta^s + \omega^l \phi^l \zeta^l ) . \end{align}

To ensure entropy production we must ensure that $\varSigma ^s + \varSigma ^l \geq 0.$ We decompose the entropy production term $\varSigma ^\varepsilon$ into contributions due to viscous heating ($\varSigma ^\varepsilon _\sigma$), internal diffusive heat flux ($\varSigma ^\varepsilon _{{k}}$) and interactions between the phases ($\varSigma ^\varepsilon _{X}$):

(2.52)\begin{equation} \varSigma^\varepsilon = \varSigma^\varepsilon_\sigma + \varSigma^\varepsilon_{{k}} + \varSigma^\varepsilon_{X}. \end{equation}

We assume that dissipative phenomena associated with only one phase depend exclusively on the variables of this phase. Thus, viscous dissipation of one phase depends only on the internal stresses in that phase and, similarly, thermal dissipation of one phase depends only on that phase's internal heat flux:

(2.53a,b)\begin{equation} \varSigma^\varepsilon_\sigma = \frac{\phi^\varepsilon \boldsymbol{\sigma}^\varepsilon : \boldsymbol{\nabla} \boldsymbol{u}^\varepsilon}{T^\varepsilon} , \quad \varSigma^\varepsilon_{{k}} =- \frac{\boldsymbol{ k }^\varepsilon \boldsymbol{\cdot} \boldsymbol{\nabla} T^\varepsilon}{T^\varepsilon} . \end{equation}

We choose the usual Newtonian form for the deviatoric stress tensor $\boldsymbol {\sigma }^\varepsilon$ and a simple description for the heat flux vector $\boldsymbol { k }^\varepsilon$:

(2.54a,b)\begin{equation} \sigma_{ij}^\varepsilon = \mu^\varepsilon \left( \frac{\partial u_i^\varepsilon}{\partial x_j} + \frac{\partial u_j^\varepsilon}{\partial x_i} - \frac{2}{3} \frac{\partial u_k^\varepsilon}{\partial x_k} \delta_{ij} \right), \quad \boldsymbol{ k }^\varepsilon =-\frac{\phi^\varepsilon k^\varepsilon}{T^\varepsilon} \boldsymbol{\nabla} T^\varepsilon , \end{equation}

where $\mu ^\varepsilon$ is the dynamic viscosity of phase $\varepsilon$ and $k^\varepsilon$ is the thermal conductivity of phase $\varepsilon$. With forms (2.54a,b), it can be easily demonstrated that both the viscous heating term and the diffusive heating term (2.53a,b) are positive definite: $\varSigma ^\varepsilon _\sigma = {(1/2) \phi ^\varepsilon \boldsymbol {\sigma }^\varepsilon : \boldsymbol {\sigma }^\varepsilon }/ T^\varepsilon \geq 0, \varSigma ^\varepsilon _{k} = \phi ^\varepsilon k^\varepsilon | (\boldsymbol {\nabla } T^\varepsilon )/{T^\varepsilon } |^2\geq 0$. Thus, to ensure entropy production, we only need to ensure that $\varSigma _{X}^s + \varSigma _{X}^l \geq 0$. It is useful to decompose the interaction term $\varSigma ^\varepsilon _{X}$ into a sum of the mean mixture entropy production $\bar {X}$ and interphase entropy exchange $\chi ^\varepsilon$ as follows:

(2.55a–c)\begin{equation} \varSigma^\varepsilon_{X} = \bar{X} + \chi^\varepsilon , \quad \bar{X} =\tfrac{1}{2} ( \varSigma^s_{X} + \varSigma^l_{X} ), \quad \chi^s =- \chi^l. \end{equation}

Using (2.55a–c), we can write $T^s \varSigma ^s_{X} + T^l \varSigma ^l_{X} = ( T^s +T^l ) \bar {X} + \Delta T \chi ^s$, and thus (2.51) becomes

(2.56)\begin{equation} \bar{X} = \frac{1}{T^s +T^l} ( - \boldsymbol{D}^s \boldsymbol{\cdot} \Delta \boldsymbol{u} - \chi^s \Delta T - \varGamma^s \Delta G_* + \Delta p ( -\omega^s \phi^s \zeta^s + \omega^l \phi^l \zeta^l ) ) . \end{equation}

Terms $\lbrace \boldsymbol {D}^s, \chi ^s, \varGamma ^s, \Delta p \rbrace$ are still unknown and require closure. At this stage we follow the standard procedure of irreversible thermodynamics (e.g. Roberts & Loper Reference Roberts and Loper1987; Ricard Reference Ricard2007; De Groot & Mazur Reference De Groot and Mazur2013) whereby the closure terms $\lbrace \boldsymbol {D}^s, \chi ^s, \varGamma ^s, \Delta p \rbrace$ are identified as thermodynamic fluxes which can be represented as linear functions of thermodynamic scalar forces. In (2.56) each thermodynamic force appears multiplied by its conjugate flux. We adopt the simplest closure approach whereby each force is linearly related to its conjugate flux:

(2.57)$$\begin{gather} \boldsymbol{D}^s =- \varLambda_D \Delta \boldsymbol{u} =- \boldsymbol{D}^l, \end{gather}$$

(2.58)$$\begin{gather}\chi^s =-\varLambda_\chi \Delta T =-\chi^l, \end{gather}$$

(2.59)$$\begin{gather}\varGamma^s =- \varLambda_\varGamma \Delta G_* =-\varGamma^l, \end{gather}$$

(2.60)$$\begin{gather}\Delta p = \varLambda_{\Delta p} ( -\omega^s \phi^s \zeta^s + \omega^l \phi^l \zeta^l ), \end{gather}$$

where $\varLambda _D, \varLambda _\chi, \varLambda _\varGamma, \varLambda _{\Delta p} > 0$ are phenomenological coefficients. Equation (2.56) can thus be written as

(2.61)\begin{equation} \bar{X} = \frac{1}{T^s +T^l} \left( \frac{| \boldsymbol{D}^s |^2}{\varLambda_D} + \frac{( \chi^s )^2}{\varLambda_\chi} + \frac{( \varGamma^s )^2}{\varLambda_\varGamma} + \frac{(\Delta p)^2}{\varLambda_{\Delta p}} \right). \end{equation}

The term $\boldsymbol {D}^s$, which has dimensions of momentum density over time, represents momentum transfer due to friction between the phases. Coefficient $\varLambda _D$, which has dimensions of density over time, is sometimes referred to as the hydraulic conductivity (Bear Reference Bear1988), and acts to reduce the velocity differences between phases. In the tight-coupling limit $\varLambda _D \to \infty$ (Roberts & Loper Reference Roberts and Loper1987), the two phases move with the same velocity ($\Delta \boldsymbol {u} = 0$). In two-phase studies of predominantly solid systems this is represented by Darcy drag or Stokes drag (see e.g. Bercovici et al. Reference Bercovici, Ricard and Schubert2001; Michaut et al. Reference Michaut, Ricard, Bercovici, Sparks and Steve2013). Following such an approach, however, results in equations which are not materially invariant. Various simple symmetric forms of the interaction coefficient $\varLambda _D$ have been used in the literature (see e.g. Bercovici & Michaut Reference Bercovici and Michaut2010; Degruyter et al. Reference Degruyter, Bachmann, Burgisser and Manga2012; Michaut et al. Reference Michaut, Ricard, Bercovici, Sparks and Steve2013). We assume that $\varLambda _D$ has the simple symmetric form

(2.62)\begin{equation} \varLambda_D = c_D \frac{\phi^s \rho^s \phi^l \rho^l}{\bar{\rho}} , \end{equation}

where $c_D$ is a constant coefficient of interphase friction (rate of frictional momentum exchange). Scaling $\varLambda _D$ with $\phi ^s\phi ^l$ ensures that momentum transfers cease in the pure-phase limit ($\phi ^\varepsilon \to 1$) (Keller & Suckale Reference Keller and Suckale2019). Note that with this parametrisation, the frictional momentum exchange becomes proportional to the phase separation flux: $\boldsymbol {D}^s = -c_D \boldsymbol {j}$ (recall (2.14)).

The entropy exchange term $\chi ^s$ is a macroscopic parametrisation of the microscale heat transfer between solid and liquid phases. Coefficient $\varLambda _\chi$ in (2.58) represents the rate of thermal equilibration between the phases. If thermal equilibration occurs via thermal conduction between solid and liquid, then equilibration of $T^s$ and $T^l$ occurs on a time scale $t_T = \bar {r}^2/\bar {\kappa }$, where $\bar {r}$ is the average solid grain radius and $\bar {\kappa }$ is a mean thermal diffusivity (Roberts & Loper Reference Roberts and Loper1987). Using the maximum grain size of 3 cm as estimated by Walker, Davies & Wilson (Reference Walker, Davies and Wilson2021), and thermal diffusivity of $10^{-5} \ {\rm m}^2 \,{\rm s}^{-1}$, the thermal equilibration time is of the order of $10\ {\rm s}$. Thus thermal equilibration between phases is possibly nearly instantaneous (at least on geophysically relevant time scales).

Equation (2.59) controls the kinetics of melting and freezing and by consequence defines the equilibrium condition. According to (2.59) freezing occurs whenever the effective Gibbs free energy of the solid is less than that of the liquid (i.e. when the solid phase is energetically favoured over the liquid). The coefficient $\varLambda _\varGamma$ in (2.59) represents the rate at which the system tends to phase equilibrium. For a pure slurry, the relaxation to phase equilibrium occurs on the same time scale as the relaxation to equality of temperatures (Roberts & Loper Reference Roberts and Loper1987).

Equation (2.60) states that the difference between mechanical pressures $\Delta p$ of two phases is a result of compaction and dilation of the phases ($\zeta ^\varepsilon$). This compaction pressure $\Delta p$ appears in the momentum equations as an additional isotropic stress (see (C3c) and (C3d)). Thus, deformation in a two-phase medium depends on two viscosities: shear viscosity $\mu ^\varepsilon$ which relates deviatoric stress and deviatoric strain rate (2.54a) (just as for a pure-phase material) and compaction viscosity $\omega ^\varepsilon \varLambda _{\Delta p}$ which describes resistance to isotropic compaction (Katz et al. Reference Katz, Jones, Rudge and Keller2022).

Coefficient $\omega ^\varepsilon, 0 \leq \omega ^\varepsilon \leq 1$, controls the partitioning of the pressure jump between the two phases and represents the fraction of the volume-averaged surface force exerted on each phase (Bercovici & Ricard Reference Bercovici and Ricard2003). At a microscopic level, the solid–liquid interfaces are not sharp discontinuities but correspond to layers (called ‘selvedge’ layers) of disorganised atom distributions. The exact value of $\omega ^\varepsilon$ is related to the microscopic behaviour of the two phases (molecular bond strengths and thickness of the interfacial selvage layers) and measures the extent to which the microscopic interface layer is embedded in one phase more than the other (Ricard Reference Ricard2007). Bercovici & Ricard (Reference Bercovici and Ricard2003) provide an expression for $\omega ^\varepsilon$ in terms of the dynamic viscosities of the phases:

(2.63)\begin{equation} \omega^\varepsilon = \frac{\phi^\varepsilon \mu^\varepsilon}{\phi^s \mu^s + \phi^l \mu^l}. \end{equation}

The exact form of the relation for $\omega ^\varepsilon$ is probably not important for many geological applications for which $\mu ^s \gg \mu ^l$. In these cases, regardless of the relation for $\omega ^\varepsilon$, one can assume $\omega ^s = 1$ and $\omega ^l = 0$. In this case, (2.60) reduces to $\Delta p = -\varLambda _{\Delta p} \phi ^s \zeta ^s$, or, without loss of generality, $\Delta p = -\mathcal {C} \zeta ^s$, where $\mathcal {C}$ is the compaction viscosity of the solid phase, which in general varies with the liquid volume fraction ($\mathcal {C}=\mathcal {C}(\phi ^l)$). It is unknown what is the appropriate form of compaction viscosity in the context of two-phase iron mixture at F-layer conditions. Roberts & Loper (Reference Roberts and Loper1987) suggest that a slurry should be in mechanical pressure equilibrium ($\Delta p = \mathcal {C}=0$). On the other hand, it is well known that in the context of partially molten rock, the compaction viscosity can be similar in magnitude to the shear viscosity. For simplicity, we first perform detailed analysis for the case with uniform shear viscosity $\mu ^\varepsilon$ (with $\mathcal {C} = 0$), and return to discuss the influence of porosity-dependent compaction viscosity in § 5.6.

The full system of governing equations consists of two equations of mass continuity (2.12), two equations of momentum (2.26a), (2.26b), two equations of entropy (2.30) and two thermodynamic relations each for density (2.45), entropy (2.46) and Gibbs free energy (2.44) (or effective Gibbs free energy (2.48)). This adds up to 16 equations for 16 variables: $\phi ^s, P, \boldsymbol {u}^s, \boldsymbol {u}^l, T^s, T^l, S^s, S^l, \rho ^s, \rho ^l, G^s, G^l$ (or $G_*^s, G_*^l$) (the full set of equations is summarised in § C.1).

2.7. Fast-melting limit

A significant simplification of the governing equations can be achieved by assuming that rates of thermal and phase equilibration are infinite: $\varLambda _\chi \to \infty, \varLambda _\varGamma \to \infty$ (Roberts & Loper Reference Roberts and Loper1987). With $\varLambda _\chi \to \infty$, the temperature of the two phases becomes equal:

(2.64)\begin{equation} T^s = T^l = T.\end{equation}

With infinitely fast melting $\varLambda _\varGamma \to \infty$, (2.59) necessitates that

(2.65)\begin{equation} \Delta G_* = 0 .\end{equation}

The equilibrium condition (2.65) must hold for all thermodynamic states; thus the total derivative of (2.65) must be zero (i.e. $\mathrm {d} G_*^s -\mathrm {d} G_*^l = 0$). Recalling (2.48a) and (2.48b), the equilibrium condition results in the following relation:

(2.66)\begin{equation} \left( \frac{1}{\rho^s} - \frac{1}{\rho^l} \right) \mathrm{d} P = \Delta S \, \mathrm{d} T - \Delta p \left( \frac{\omega^s}{( \rho^s )^2}\, \mathrm{d} \rho^s + \frac{\omega^l}{( \rho^l )^2} \,\mathrm{d} \rho^l \right).\end{equation}

Constraint (2.66) means that the temperature $T$ and pressure $P$ are no longer independent variables. The thermodynamic equilibrium condition (2.65) does not mean that phase change is not allowed. Even at equilibrium $\varGamma ^s$ can be finite, though it may no longer be determined directly through (2.59). Instead, $\varGamma ^s$ is determined indirectly by the constraint (2.66) which restricts the temperature to the melting/freezing temperature. Note that in the absence of compaction ($\Delta p =0$) the fast-melting constraint (2.66) becomes identical to the classic Clapeyron relation. When the compaction pressure is non-zero ($\Delta p \neq 0$), the pressure difference between phases affects the melting temperature if the phases are compressible ($\mathrm {d} \rho ^s \neq 0, \mathrm {d} \rho ^l \neq 0$).

When the two phases share the same temperature there is no need for two separate entropy equations. Instead we may formulate an evolution equation for entropy density of the mixture $\overline {\rho S}$ (2.9). Note first that the mixture entropy flux $\overline {\rho S\boldsymbol {u}} = (\overline {\rho S})( \overline {\rho \boldsymbol {u}})/\bar {\rho }$ can be written as

(2.67)\begin{equation} \overline{\rho S\boldsymbol{u}} = \phi^s\rho^s S^s \boldsymbol{u}^s +\phi^l\rho^l S^l \boldsymbol{u}^l - \Delta S \boldsymbol{j} .\end{equation}

Summing entropy equations of individual phases (2.30), and using (2.67), we obtain the conservation equation for the entropy density of the mixture $\overline {\rho S}$,

(2.68)\begin{equation} \frac{\partial }{\partial t} ( \overline{\rho S} ) + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \overline{\rho S\boldsymbol{u}} ) + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \Delta S \boldsymbol{j} ) + \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{ k } = \varSigma,\end{equation}

with total heat flux vector $\boldsymbol { k }$ and total entropy production term $\varSigma$:

(2.69)$$\begin{gather} \boldsymbol{ k } = \boldsymbol{ k }^s + \boldsymbol{ k }^l =-\frac{\bar{k}}{T} \boldsymbol{\nabla} T, \quad \text{where} \ \bar{k} = {\phi^s k^s + \phi^l k^l}, \end{gather}$$

(2.70)$$\begin{gather}\varSigma = \varSigma^s + \varSigma^l = \frac{1}{T} \left( - \boldsymbol{ k } \boldsymbol{\cdot} \boldsymbol{\nabla} T + \phi^s \boldsymbol{\sigma}^s : \boldsymbol{\nabla} {\boldsymbol{u}^s} + \phi^l \boldsymbol{\sigma}^l : \boldsymbol{\nabla} {\boldsymbol{u}^l} +\varLambda_D | \Delta \boldsymbol{u} |^2 + \frac{( \Delta p )^2}{\varLambda_{\Delta p}} \right) . \end{gather}$$

Using the total mass conservation equation (2.13), and the definitions of $\boldsymbol { k }$ and $\varSigma$, the barycentric entropy equation (2.68) becomes

(2.71)\begin{equation} \bar{\rho} T \frac{\bar{\mathrm{D}} \bar{S}}{\bar{\mathrm{D}} t} - T \boldsymbol{\nabla}\boldsymbol{\cdot} \left( \frac{L }{T } \boldsymbol{j} \right) = \boldsymbol{\nabla}\boldsymbol{\cdot} ( \bar{k} \boldsymbol{\nabla} T ) + \varPsi,\end{equation}

where $L$ is latent heat of fusion and $\varPsi$ is deformational work:

(2.72)$$\begin{gather} L \equiv-T \Delta S, \end{gather}$$

(2.73)$$\begin{gather}\varPsi = \phi^s \boldsymbol{\sigma}^s : \boldsymbol{\nabla} {\boldsymbol{u}^s} + \phi^l \boldsymbol{\sigma}^l : \boldsymbol{\nabla} {\boldsymbol{u}^l} +\varLambda_D | \Delta \boldsymbol{u} |^2 + \frac{( \Delta p )^2}{\varLambda_{\Delta p}} . \end{gather}$$

Throughout, we shall assume that $L$ is constant.

We may convert the equation for mixture entropy $\bar {S}$ into an equation for the temperature $T$. To do that, we first need to derive the relation for the thermodynamic derivative of mixture entropy. Recall that, by definition (2.9) (and (2.6)), $\bar {S} = \psi ^s S^s + \psi ^l S^l$, which, on differentiation, gives

(2.74)\begin{equation} \mathrm{d} \bar{S} = \psi^s \mathrm{d} S^s + (1-\psi^s) \mathrm{d} S^l + \Delta S \mathrm{d} \psi^s. \end{equation}

Substituting the definition of $\mathrm {d} S^\varepsilon$ (2.46) into (2.74) we obtain the differential of mixture entropy,

(2.75)\begin{equation} \mathrm{d} \bar{S} =-\frac{\bar{\alpha}}{\bar{\rho}} \mathrm{d} P + \frac{\bar{c}_p}{T} \mathrm{d} T - \frac{L}{T} \mathrm{d} \psi^s,\end{equation}

where $\bar {\alpha }$ is the thermal expansion coefficient of the mixture, $\bar {c}_p$ is the specific heat capacity of the mixture and $L$ is latent heat of fusion per unit mass:

(2.76)$$\begin{gather} \bar{\alpha} = \phi^s \alpha^s + \phi^l \alpha^l, \end{gather}$$

(2.77)$$\begin{gather}\bar{c}_p = ( \phi^s \rho^s c_p^s + \phi^l \rho^l c_p^l )/\bar{\rho} . \end{gather}$$

Using (2.75) we can write

(2.78)\begin{equation} \bar{\rho} T \frac{\bar{\mathrm{D}} \bar{S}}{\bar{\mathrm{D}} t} =-\bar{\alpha} T \frac{\bar{\mathrm{D}} P}{\bar{\mathrm{D}} t} + \bar{c}_p \bar{\rho} \frac{\bar{\mathrm{D}} T}{\bar{\mathrm{D}} t} - L \bar{\rho} \frac{\bar{\mathrm{D}} \psi^s}{\bar{\mathrm{D}} t}.\end{equation}

On using (2.78), together with the equation for ${\bar {\mathrm {D}} \psi ^s}/{\bar {\mathrm {D}} t}$ (2.14), and the definition of latent heat (2.72), the entropy equation (2.71) becomes the temperature equation:

(2.79)\begin{equation} \bar{c}_p \bar{\rho} \frac{\bar{\mathrm{D}} T}{\bar{\mathrm{D}} t} -\bar{\alpha} T \frac{\bar{\mathrm{D}} P}{\bar{\mathrm{D}} t} -L \varGamma^s + \frac{L }{ T} \boldsymbol{j} \boldsymbol{\cdot} \boldsymbol{\nabla} T = \boldsymbol{\nabla}\boldsymbol{\cdot} ( \bar{k} \boldsymbol{\nabla} T ) + \varPsi. \end{equation}

The full system of equations governing the evolution of the fast-melting slurry consists of two equations of mass continuity (2.12), two equations of momentum (2.26a), (2.26b), the equation for the temperature of the mixture (2.79), the liquidus constraint (2.66) and two thermodynamic relations for solid and liquid densities (2.45). This adds up to 12 equations for 12 variables: $\phi ^s, P, \boldsymbol {u}^s, \boldsymbol {u}^l, T, \rho ^s, \rho ^l, \varGamma ^s$ (the full set of equations of fast-melting slurry is summarised in § C.2).

Notably, variables $S^s, S^l, G^s, G^l$, and thus (2.46), (2.44), are no longer relevant in the fast-melting limit. Similarly, $\varGamma ^s$ is no longer determined directly by (2.59). Instead we may think that the temperature $T$ follows the liquidus (2.66), and the freezing rate $\varGamma ^s$ is determined by the heat equation (2.79).

2.8. Boundary conditions

To complete the theory it is necessary to specify boundary conditions on bounding surfaces. In general the relevant conditions at both the top and bottom of the slurry are not known because they are not constrained by observations. There are a limited number of conditions that can be deduced from the governing equations, which need to be supplemented with additional assumptions and approximations based on physical intuition or experimental insight where available.

In the most general case the position of bounding surfaces may be not fixed and may evolve in time. For example, in Earth's core, over time the ICB advances in the radial direction due to inner core freezing. Thus, a boundary $B$ may move at an unknown velocity $\boldsymbol {U}_B$, which does not in general coincide with the flow velocity $\boldsymbol {u}^\varepsilon$. By integrating the governing equations over a pillbox volume that straddles the boundary and moves together with it we can derive continuity conditions that must hold at the interface. Applying this procedure to equations of mass (2.12), momentum (2.16) and total energy (2.37) leads to the following conditions:

(2.80)$$\begin{gather}{}[\![ { \phi^\varepsilon \rho^\varepsilon( \boldsymbol{u}^\varepsilon -\boldsymbol{U}_B ) \boldsymbol{\cdot} \hat{\boldsymbol{ n }} } ]\!] = 0 , \quad \text{on} \ B, \end{gather}$$

(2.81)$$\begin{gather}{}[\![ { \phi^\varepsilon \rho^\varepsilon ( \boldsymbol{u}^\varepsilon\boldsymbol{\cdot} \hat{\boldsymbol{ n }} ) ( ( \boldsymbol{u}^\varepsilon - \boldsymbol{U}_B ) \boldsymbol{\cdot} \hat{\boldsymbol{ n }} ) + \phi^\varepsilon p^\varepsilon - \phi^\varepsilon\hat{\boldsymbol{ n }}\boldsymbol{\cdot}\boldsymbol{\sigma}^\varepsilon \boldsymbol{\cdot} \hat{\boldsymbol{ n }} } ]\!] = 0, \quad \text{on}\ B, \end{gather}$$

(2.82)$$\begin{gather}{}[\![ { \phi^\varepsilon \rho^\varepsilon ( \boldsymbol{u}^\varepsilon\boldsymbol{\cdot} \hat{\boldsymbol{ t }} ) ( ( \boldsymbol{u}^\varepsilon - \boldsymbol{U}_B ) \boldsymbol{\cdot} \hat{\boldsymbol{ n }} ) - \phi^\varepsilon\hat{\boldsymbol{ t }}\boldsymbol{\cdot} \boldsymbol{\sigma}^\varepsilon \boldsymbol{\cdot} \hat{\boldsymbol{ n }} } ]\!] = 0, \quad \text{on } B, \end{gather}$$

(2.83)$$\begin{gather}{}[\![ { \overline{ \rho \mathcal{E}} ( \boldsymbol{\bar{u}} - \boldsymbol{U}_B ) + \Delta \mathcal{E} \boldsymbol{j} + \phi^s p^s \boldsymbol{u}^s + \phi^l p^l \boldsymbol{u}^l - \phi^s \boldsymbol{u}^s \boldsymbol{\cdot} \boldsymbol{\sigma}^s - \phi^l \boldsymbol{u}^l \boldsymbol{\cdot} \boldsymbol{\sigma}^l + T \boldsymbol{ \bar{k} }} ]\!] \boldsymbol{\cdot} \hat{\boldsymbol{ n }} = 0, \quad \text{on}\ B, \end{gather}$$

where $\hat {\boldsymbol { n }}$ ($\hat {\boldsymbol { t }}$) is a unit vector normal (tangent) to the boundary surface $B$ and $[\![ { A } ]\!]$ denotes the jump in a quantity $A$ across the boundary $B$ (detailed derivation of these conditions is included in Appendix D). Conditions (2.80)–(2.83) represent continuity of mass flux, tangential and normal components of momentum flux and total energy flux, respectively.

Condition (2.80) is known as the kinematic condition and specifies the relationship between normal components of velocity. Following Leal (Reference Leal2007) and Loper & Roberts (Reference Loper and Roberts1978) we additionally impose the standard dynamic condition of continuity of tangential component of velocity,

(2.84)\begin{equation} [\![ \boldsymbol{u}^\varepsilon \boldsymbol{\cdot} \hat{\boldsymbol{ t }} ]\!] = 0 , \quad \text{on} \ B.\end{equation}

This is usually known as the no-slip condition, and it is generally accepted to apply under almost all circumstances for Newtonian fluids either at solid surfaces or at fluid–fluid interfaces (Leal Reference Leal2007). Note that with (2.84) and (2.80), continuity of tangential momentum flux (2.82) reduces to the condition of continuity of tangential components of deviatoric stress, i.e. $[\![ \phi ^\varepsilon \hat {\boldsymbol { t }}\boldsymbol {\cdot } \boldsymbol {\sigma }^\varepsilon \boldsymbol {\cdot } \hat {\boldsymbol { n }} ]\!] = 0$ on $B$.

In addition to the above conditions, it is often assumed that the bounding surface is maintained at thermal equilibrium (e.g. Leal Reference Leal2007), and so the temperature varies continuously across the interface:

(2.85)\begin{equation} [\![ T ]\!] = 0, \quad \text{on} \ B.\end{equation}

Similarly, it is often assumed that the (thermodynamic) pressure is continuous on the interface (Loper & Roberts Reference Loper and Roberts1978):

(2.86)\begin{equation} [\![ P ]\!] = 0 , \quad \text{on} \ B.\end{equation}

In light of continuity of temperature (2.85) and pressure (2.86), it follows from the equation of state (2.45) that the specific density of either phase will also be continuous:

(2.87)\begin{equation} [\![ \rho^\varepsilon ]\!] = 0 , \quad \text{on} \ B.\end{equation}

Note that this does not mean that the total (mixture) density $\bar {\rho }$ is continuous across the interface since, in general, $\phi ^s$ is not necessarily continuous.

Equations (2.80)–(2.87) are continuity conditions, and thus require knowledge about variables on either side of the domain boundary. Since the values that these variables take on the exterior of the boundaries are not determined by the interior solution, additional assumptions about properties of the bounding regions are required to turn the continuity conditions into boundary conditions. Specific boundary conditions are developed presently (in § 3.1). Note also that, in general, $\boldsymbol {U}_B$ is unknown and an additional equation is required to determine it.

3.1. Problem formulation

We work in Cartesian geometry, in a plane parallel layer of depth $d$ (figure 1). Coordinates $x, y$ denote horizontal directions and $z$ denotes the vertical. Gravity acts in the vertical direction $\boldsymbol { g }=-g\hat {\boldsymbol { e }}_z$. The bottom boundary ($z=0$) is representative of the ICB, while the top boundary ($z=d$) marks the end of the F-layer. For simplicity, we assume that the boundaries are fixed and not moving ($\boldsymbol {U}_B = 0$).

Figure 1. Schematic of the problem domain.

The problem of a single-phase compressible fluid bounded in the vertical direction requires nine boundary conditions: one on the density, two on the temperature and two conditions (for each bounded dimension) on each of the three components of the velocity vector. In the two-phase problem, we have an additional density field and an additional velocity vector, and thus we have to include seven additional boundary conditions, which leads to a total of 16 boundary conditions. In practice, there are more continuity conditions than the number of boundary conditions that can be applied; only a subset of these conditions can be used to develop practicable boundary conditions. The choice of which continuity conditions should take priority depends on the nature of the bounding region.

We set the values of density and temperature at the bottom boundary to their reference values (2.85), (2.87):

(3.1a–c)\begin{equation} \rho^l = \rho_r^l, \quad \rho^s = \rho_r^s, \quad T = T_r \quad \text{on} \ z = 0.\end{equation}

We assume that the inner core is static and impermeable. Conditions of continuity of mass flow (2.80) and tangential velocity (2.84) then reduce to

(3.2)\begin{equation} \boldsymbol{u}^s = \boldsymbol{u}^l = 0 \quad \text{on} \ z = 0, \end{equation}

while the condition of continuity of energy flux (2.83) reduces to the statement of continuity of heat flux:

(3.3)\begin{equation} \bar{k}\frac{\partial T}{\partial z} =- q \quad \text{on} \ z = 0, \end{equation}

where $q$ is the secular cooling (per unit area) of the inner core and $\bar {k}$ is conductivity of the slurry (recall (2.69)).

Unlike the rigid surface on the bottom, the top boundary represents the interface between the F-layer slurry and the rest of the outer core. In the F-layer, the temperature is constrained to the liquidus and the solid volume fraction is non-zero everywhere (including at the top boundary). On the other hand, the convection zone is entirely liquid and the temperature does not follow the liquidus. To reconcile these two regions, we envisage a thin transition region between the F-layer and convection zone proper (see figure 1), where the temperature profile departs from the liquidus and the solid fraction diminishes to zero. Thus we assume that $\phi ^\varepsilon$ is continuous across $z = d$. In that case, the conditions of continuity of tangential and normal stress (2.81), (2.82) reduce to

(3.4a–c)\begin{equation} [\![ { \sigma^\varepsilon_{xz} } ]\!]= 0, \quad [\![ { \sigma^\varepsilon_{yz} } ]\!]= 0, \quad [\![ { \sigma^\varepsilon_{zz} } ]\!]= 0 \quad \text{on} \ z = d. \end{equation}

Since we do not solve for the velocities beyond $z=d$ we cannot directly use the above conditions. Instead, we make the simplifying assumption that top boundary is a stress-free surface where all components of the deviatoric stress tensor are zero:

(3.5a–c)\begin{equation} \sigma^\varepsilon_{xz} = 0, \quad \sigma^\varepsilon_{yz} = 0, \quad \sigma^\varepsilon_{zz} = 0 \quad \text{on} \ z = d.\end{equation}

Equations (3.1a–c)–(3.3) and (3.5a–c) form a complete set of 16 boundary conditions necessary for the solution of the two-phase slurry problem.

3.2. Dimensionless equations

We begin by non-dimensionalising the governing equations (2.12), (2.26a), (2.26b), (2.79), (2.66), (2.45) (or equivalently, (C3a)–(C3e)). For units of density (liquid and solid) and temperature we choose their values at the bottom of the layer ($z=0$): $\rho _r^l, \rho _r^s$ and $T_r$, respectively. (The subscript ‘$r$’ is used throughout to denote representative/reference values of quantities.) We scale length with the depth of the layer $d$, time with the gravitational free-fall time scale $\sqrt {d/g}$, velocity with $\sqrt {g d}$, pressure with $\rho _r^l gd$ and the mass production term $\varGamma ^s$ with $\rho _r^l \sqrt {g / d}$:

(3.6)\begin{equation} \left.\begin{array}{c@{}} \boldsymbol{ x } = {\rm d} \hat{\boldsymbol{ x }}, \quad t = ( {d}/{g} )^{1/2} \hat{t}, \quad \boldsymbol{u}^\varepsilon = ( g d )^{1/2} \hat{\boldsymbol{u}}^\varepsilon, \quad P = ( \rho_r^l gd ) \hat{P} ,\\ \rho^l = \rho_r^l \hat{\rho}^l, \quad \rho^s = \rho_r^s \hat{\rho}^s, \quad T = T_r \hat{T}, \quad \varGamma^s = \rho_r^l ( g/d )^{1/2} \hat{\varGamma}^s . \end{array}\right\} \end{equation}

(Note that volume fractions $\phi ^\varepsilon$ are dimensionless and bounded between 0 and 1, and thus do not require any scaling.) We assume that material coefficients of individual phases $\alpha ^s, \alpha ^l, \beta ^s, \beta ^l, c_p^s, c_p^l, k^s, k^l, \mu ^s, \mu ^l$ are constant. We also set $\omega ^s =1, \omega ^l = 0$ and $\gamma ^s = \gamma ^l = 1/2$. Applying outlined scalings (3.6) to the governing equations, and omitting hats for ease of notation, we obtain the following set of dimensionless equations:

(3.7a)$$\begin{gather} \lambda_\rho \left( \frac{\partial (\phi^s \rho^s)}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^s \rho^s \boldsymbol{u}^s ) \right) = \varGamma^s, \end{gather}$$

(3.7b)$$\begin{gather}\frac{\partial (\phi^l \rho^l)}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^l \rho^l \boldsymbol{u}^l ) =-\varGamma^s, \end{gather}$$

(3.7c)\begin{gather} \lambda_\rho \phi^s \rho^s \left( \frac{\partial \boldsymbol{u}^s}{\partial t} + \boldsymbol{u}^s \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}^s \right) =- \phi^s \boldsymbol{\nabla} P - \lambda_\rho \phi^s \rho^s \boldsymbol{\hat{e}}_z - K \boldsymbol{j} - \frac{1}{2} \varGamma^s \Delta \boldsymbol{u} \nonumber\\ \quad + \lambda_\mu \sqrt{ \frac{Pr}{R}} \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^s ( \boldsymbol{\sigma}^s + \mathcal{C} (\phi^l) \zeta^s \boldsymbol{ I } ) ) , \end{gather}

(3.7d)\begin{gather} \phi^l \rho^l \left( \frac{\partial \boldsymbol{u}^l}{\partial t} + \boldsymbol{u}^l \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}^l \right) =- \phi^l \boldsymbol{\nabla} P - \phi^l \rho^l \boldsymbol{\hat{e}}_z + K \boldsymbol{j} - \frac{1}{2} \varGamma^s \Delta \boldsymbol{u} + \sqrt{ \frac{Pr}{R}} \boldsymbol{\nabla}\boldsymbol{\cdot} ( \phi^l \boldsymbol{ \sigma }^l ) , \end{gather}

(3.7e)\begin{gather} \bar{c}_p \bar{\rho} \left( \frac{\partial T}{\partial t} + \boldsymbol{\bar{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} T \right) - \alpha^* D \bar{\alpha} T \left( \frac{\partial P}{\partial t} + \boldsymbol{\bar{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} P \right) - \frac{\varGamma^s}{St} + \frac{\lambda_\rho}{St} \frac{ \boldsymbol{j} \boldsymbol{\cdot} \boldsymbol{\nabla} T}{ T} = \frac{1}{\sqrt{R Pr}} \boldsymbol{\nabla}\boldsymbol{\cdot} ( \bar{k} \boldsymbol{\nabla} T )\nonumber\\ \quad + D K \boldsymbol{j} \boldsymbol{\cdot} \Delta \boldsymbol{u} +D \sqrt{\frac{Pr}{R} } \left( \frac{1}{2} \lambda_\mu \phi^s \boldsymbol{\sigma}^s : \boldsymbol{\sigma}^s + \frac{1}{2} \phi^l \boldsymbol{\sigma}^l : \boldsymbol{\sigma}^l + \lambda_\mu \phi^s \mathcal{C}(\phi^l) ( {\zeta}^s )^2 \right) , \end{gather}

(3.7f)$$\begin{gather} \left( \frac{1}{\rho^l} - \frac{1}{\lambda_\rho \rho^s} \right) \mathrm{d} P =\frac{1}{St\,D}\frac{1}{T} \, \mathrm{d} T - \frac{ \lambda_\mu }{\lambda_\rho} \sqrt{\frac{Pr}{R} } \frac{\mathcal{C}(\phi^l) \zeta^s}{ ( \rho^s )^2} \,\mathrm{d} \rho^s , \end{gather}$$

(3.7g)$$\begin{gather}\mathrm{d} \rho^s = \rho^s ( \lambda_\beta \beta^* \mathrm{d} P - \lambda_\alpha \alpha^* \mathrm{d} T ), \end{gather}$$

(3.7h)$$\begin{gather}\mathrm{d} \rho^l = \rho^l ( \beta^* \mathrm{d} P - \alpha^* \mathrm{d} T ), \end{gather}$$

where

(3.8) \begin{equation} \left.\begin{array}{c@{}} \displaystyle \phi^l = 1 - \phi^s, \quad \Delta \boldsymbol{u} = \boldsymbol{u}^s-\boldsymbol{u}^l, \quad \boldsymbol{j} = \dfrac{\phi^s \rho^s \phi^l \rho^l}{\bar{\rho}}\Delta \boldsymbol{u}, \\ \displaystyle \bar{\rho} = \lambda_\rho \phi^s \rho^s + \phi^l \rho^l , \quad \bar{c}_p = ( \lambda_{c_p} \lambda_\rho \phi^s \rho^s + \phi^l \rho^l )/\bar{\rho} ,\quad \bar{\alpha} = \lambda_\alpha \phi^s + \phi^l ,\\ \displaystyle \bar{k} = \lambda_k \phi^s + \phi^l ,\quad \boldsymbol{\bar{u}} = ( \lambda_\rho \phi^s\rho^s \boldsymbol{u}^s + \phi^l \rho^l \boldsymbol{u}^l )/\bar{\rho} ,\\ \displaystyle \sigma_{ij}^\varepsilon = \dfrac{\partial u_i^\varepsilon}{\partial x_j} + \dfrac{\partial u_j^\varepsilon}{\partial x_i} - \dfrac{2}{3} \dfrac{\partial u_k^\varepsilon}{\partial x_k} \delta_{ij}, \quad \zeta^s = \dfrac{1}{\rho^s}\left( \dfrac{\partial \rho^s}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} (\rho^s \boldsymbol{u}^s) \right). \end{array}\right\} \end{equation}

The dimensionless numbers are

(3.9)\begin{align} \left.\begin{array}{c@{}} \displaystyle R = \dfrac{c_{p}^l ( \rho_r^l )^2 g d^3}{k^l \mu^l}, \quad Pr = \dfrac{c_{p}^l \mu^l}{k^l}, \quad D = \dfrac{gd}{c_{p}^l T_r }, \quad St = \dfrac{c_{p}^l T_r}{L}, \quad K = \dfrac{c_D \rho_r^s }{\rho_r^l (g/d)^{1/2}}, \\ \displaystyle \theta = \dfrac{q d}{k_r^l T_r}, \quad \alpha^* = \alpha^l T_r, \quad \beta^* = \beta^l \rho_r^l g d,\\ \displaystyle \lambda_\alpha = \dfrac{\alpha^s}{\alpha^l}, \quad \lambda_\beta = \dfrac{\beta^s}{\beta^l}, \quad \lambda_{c_p} = \dfrac{c_{p}^s}{c_{p}^l}, \quad \lambda_\rho = \dfrac{\rho_{r}^s}{\rho_{r}^l}, \quad \lambda_k = \dfrac{k^s}{k^l}, \quad \lambda_\mu = \dfrac{\mu^s}{\mu^l}. \end{array}\right\} \end{align}

Parameter $R$ quantifies the strength of the gravitational force as measured against viscous and thermal diffusive forces (it is similar to the Rayleigh number associated with buoyancy-driven flow although this analogy is not complete since $R$ as defined here is missing a factor describing the density difference across the layer); the Prandtl number $Pr$ is the ratio of viscosity to thermal diffusivity; the dissipation number $D$ measures the influence of compressibility/stratification; the Stefan number $St$ measures the ratio between sensible and latent heat; $K$ describes the effects of interphase drag/friction; parameter $\theta$ arises from the heat-flux boundary condition (3.2), which in dimensionless form becomes $\mathrm {d} {T}/\mathrm {d}{z} = - \theta /\bar {k}$, where $\theta = q d/k_r^l T_r$ is the dimensionless heat flux (or temperature gradient) at the bottom boundary; $\alpha ^*$ and $\beta ^*$ are non-dimensionalised coefficients of thermal expansion and isothermal compression respectively; parameters $\lbrace \lambda _\alpha, \lambda _\beta, \lambda _{c_p}, \lambda _\rho, \lambda _k, \lambda _\mu \rbrace$ measure ratios of physical properties between solid and liquid phases. We note that studies of two-phase flows in the Earth's mantle often employ an alternative non-dimensionalisation based on compaction length (given by the square root of the ratio of viscous force and interphase drag) and Darcy velocity (ratio of relative buoyancy and interphase drag) (see e.g. McKenzie Reference McKenzie1984).

Table 1 provides estimates of dimensional parameters relevant for Earth's core, while estimates of the dimensionless parameters are included in table 2. A detailed discussion of geophysically relevant ranges for many of the parameters is given in Wong et al. (Reference Wong, Davies and Jones2021). In order to reduce the parameter space that needs to be explored we keep the following parameters fixed to values that are relevant to Earth's F-layer:

(3.10)\begin{equation} \left.\begin{array}{c@{}} \alpha^* = 0.05, \quad \beta^* = 0.005, \quad D = 0.1, \quad Pr = 0.1, \quad St = 4,\\ \theta = 0.09, \quad \lambda_\rho = 1.25, \quad \lambda_\alpha = \lambda_\beta = \lambda_{c_p} = \lambda_k = 1. \end{array}\right\} \end{equation}

These parameters are reasonably well known and we expect that uncertainties on their values would have only a minor effect on the solutions presented below. It is worth noting that setting $\lambda _\alpha = \lambda _{c_p} = \lambda _k = 1$ simplifies the temperature equation (3.7e) and the heat flux boundary condition (${\mathrm {d} T}/{\mathrm {d} z} = - \theta /\bar {k}$) because the mixture thermodynamic coefficients become constant: $\bar {\alpha } = \bar {c}_p \bar {\rho } = \bar {k} = 1$. By contrast, the magnitudes of $R$ and $\lambda _{\mu }$ are numerically intractable and so we can only search for limiting behaviour as these parameters are increased towards geophysical values.

Table 1. Estimates of physical parameters characteristic of pure iron at ICB conditions.

Table 2. Estimates of dimensionless numbers at core conditions (based on table 1), and values used in this study.

The value of $K$ is hard to estimate for the F-layer slurry; it depends on the momentum exchange coefficient $c_D$ introduced through the parametrisation of the drag term $\varLambda _D = c_D \phi ^s \rho ^s \phi ^l \rho ^l/\bar {\rho }$ (2.57). In magmatic settings the parametrisation of $\varLambda _D$ is based on the permeability of the solid matrix, $\varLambda _D^{m} = \mu ^l/(k_\phi {\phi ^l}^{n-2})$, where $n=2$ is usually taken and $k_\phi$ is a permeability constant (units $\textrm {m}^2$) whose values can range from $10^{-7}$ (Spencer, Katz & Hewitt Reference Spencer, Katz and Hewitt2020) to $5\times 10^{-10}$ (Šrámek et al. Reference Šrámek, Ricard and Bercovici2007). To match the magnitude of the drag term, $\varLambda _D \sim \varLambda _D^{m}$, we may approximate $c_D = \mu ^l/( k_\phi \rho _r^l)$. Using $k_\phi = 5\times 10^{-10}$–$10^{-7}$ gives $c_D = 8$–$1600\ \textrm {s}^{-1}$, and the dimensionless interaction parameter in the range $K = 10^3$–$3 \times 10^5$. However, such low values of $k_\phi$ characterise scenarios where the liquid volume fraction is small and the solid phase forms an interlocking matrix, and are therefore not appropriate for the slurry F-layer where the solid volume fraction must be very small (of the order of 1 %; Gubbins et al. Reference Gubbins, Masters and Nimmo2008). We expect the effective permeability of the slurry to be somewhat larger, implying smaller $K$. In the absence of experimental information at planetary core conditions, we explore the range $K = 10^2$–$10^4$, which corresponds to $c_D = 0.5$–$50\ \textrm {s}^{-1}$.

4.1. Governing equations and boundary conditions

We seek a time-independent solution that depends on the $z$ coordinate only, with unidirectional motion along the vertical direction (no horizontal motion), i.e. $\boldsymbol {u}^\varepsilon = ( 0, \, 0, \, u_z^\varepsilon (z) )$. In this 1-D state the horizontal components of velocity of each phase drop out, along with equations of horizontal momenta, and the system of (3.7a)–(3.7f) reduces from 12 to 8 equations:

(4.1a)$$\begin{gather} \lambda_\rho \frac{\mathrm{d} }{\mathrm{d} z} ( \phi^s \rho^s u_z^s ) = \varGamma^s, \end{gather}$$

(4.1b)$$\begin{gather}\frac{\mathrm{d} }{\mathrm{d} z} ( \phi^l \rho^l u_z^l ) =-\varGamma^s, \end{gather}$$

(4.1c) $$\begin{gather}\lambda_\rho \phi^s \rho^s u_z^s \frac{\mathrm{d} u_z^s}{\mathrm{d} z} =- \phi^s \frac{\mathrm{d} P}{\mathrm{d} z} - \lambda_\rho \phi^s \rho^s - Kj_z - \frac{1}{2}\varGamma^s {\Delta u_z} \notag\\ + \lambda_\mu \sqrt{ \frac{Pr}{R}} \frac{\mathrm{d} }{\mathrm{d} z} \left( \frac{4}{3} \phi^s \frac{\mathrm{d} u_z^s}{\mathrm{d} z} + \phi^s \mathcal{C} (\phi^l) \zeta^s \right), \end{gather}$$

(4.1d)$$\begin{gather}\phi^l \rho^l u_z^l \frac{\mathrm{d} u_z^l}{\mathrm{d} z} =- \phi^l \frac{\mathrm{d} P}{\mathrm{d} z} - \phi^l \rho^l + Kj_z - \frac{1}{2}\varGamma^s \Delta u_z + \frac{4}{3} \sqrt{ \frac{Pr}{R}} \frac{\mathrm{d} }{\mathrm{d} z} \left( \phi^l \frac{\mathrm{d} u_z^l}{\mathrm{d} z} \right), \end{gather}$$

(4.1e)\begin{gather} \bar{c}_p \bar{\rho} \bar{u}_z \frac{\mathrm{d} T}{\mathrm{d} z} -\alpha^* D \bar{\alpha} T \bar{u}_z \frac{\mathrm{d} P}{\mathrm{d} z} - \frac{\varGamma^s}{St} + \frac{\lambda_\rho}{St} \frac{j_z }{T} \frac{\mathrm{d} T}{\mathrm{d} z} = \frac{1}{\sqrt{R\, Pr}} \frac{\mathrm{d} }{\mathrm{d} z} \left( \bar{k} \frac{\mathrm{d} T}{\mathrm{d} z} \right) + D \, K j_z \Delta u_z \nonumber\\ \quad + D \sqrt{\frac{Pr}{R}} \left( \frac{4 }{3}\lambda_\mu \phi^s \left( \frac{\mathrm{d} u_z^s}{\mathrm{d} z} \right)^2 + \frac{4 }{3} \phi^l \left( \frac{\mathrm{d} u_z^l}{\mathrm{d} z} \right)^2 + \lambda_\mu \phi^s \mathcal{C} (\phi^l) ( {\zeta}^s )^2 \right), \end{gather}

(4.1f)$$\begin{gather} \left( \frac{1}{\rho^l} - \frac{1}{\lambda_\rho \rho^s} \right) \frac{\mathrm{d} P}{\mathrm{d} z} =\frac{1}{St\,D}\frac{1}{T} \, \frac{\mathrm{d} T}{\mathrm{d} z} - \frac{ \lambda_\mu }{\lambda_\rho} \sqrt{\frac{Pr}{R}} \frac{\mathcal{C}(\phi^l) \zeta^s}{ ( \rho^s )^2} \frac{\mathrm{d} \rho^s}{\mathrm{d} z} , \end{gather}$$

(4.1g)$$\begin{gather}\frac{\mathrm{d} \rho^s }{\mathrm{d} z} = \rho^s \left( \lambda_\beta \beta^* \frac{\mathrm{d} P}{\mathrm{d} z} - \lambda_\alpha \alpha^* \frac{\mathrm{d} T}{\mathrm{d} z} \right) , \end{gather}$$

(4.1h)$$\begin{gather}\frac{\mathrm{d} \rho^l}{\mathrm{d} z} = \rho^l \left( \beta^* \frac{\mathrm{d} P}{\mathrm{d} z} - \alpha^* \frac{\mathrm{d} T}{\mathrm{d} z} \right). \end{gather}$$

The eight variables are: $\phi ^s, P, u_z^s, u_z^l, T, \rho ^s, \rho ^l, \varGamma ^s$. Auxiliary relations included in (3.8) give $\phi ^l = 1 - \phi ^s, \zeta ^s = {\mathrm {d} u_z^s}/{\mathrm {d} z} + (u_z^s/\rho ^s) {\mathrm {d} \rho ^s}/{\mathrm {d} z}, j_z = \phi ^s\rho ^s \phi ^l \rho ^l \Delta u_z /\bar {\rho }$, etc. Note also that four boundary conditions on the horizontal components of phase velocities (3.3) and four boundary conditions on tangential stresses (3.5a–c) are automatically satisfied in the 1-D state, and thus the number of boundary conditions reduces from 16 to 8. Boundary conditions (3.1a–c), (3.2), (3.3), (3.5a–c), in dimensionless form, read

(4.2)\begin{equation} \left.\begin{array}{c@{}} \displaystyle \rho^s = 1, \quad \rho^l = 1, \quad T = 1, \quad \dfrac{\mathrm{d} T}{\mathrm{d} z} =- \theta, \quad u_z^s = u_z^l = 0 \quad \text{on } z= 0;\\ \displaystyle \dfrac{\mathrm{d} u_z^s}{\mathrm{d} z} = \dfrac{\mathrm{d} u_z^l}{\mathrm{d} z} = 0 \quad \text{on } z = 1. \end{array} \right\} \end{equation}

4.2. Characteristics of the 1-D time-independent state

The 1-D time-independent solution of the two-phase system cannot be static ($u_z^s =u_z^l=0$); if that were the case the governing equations reduce to an overprescribed (contradictory) system. With $u_z^s =u_z^l=0$, solid and liquid momentum equations (4.1c), (4.1d) give

(4.3a,b)\begin{equation} 0 =- \phi^s \frac{\mathrm{d} P}{\mathrm{d} z} - \lambda_\rho \phi^s \rho^s; \quad 0 =- \phi^l \frac{\mathrm{d} P}{\mathrm{d} z} - \phi^l \rho^l . \end{equation}

Unless either $\phi ^s = 0$ or $\phi ^l=0$ (which is unlike the slurry state where solid fraction varies $0<\phi ^s<1$) these two equations cannot be satisfied if the two phases have different densities ($\lambda _\rho \rho ^s \neq \rho ^l$). Thus, in general the basic state flow velocities have to be non-zero.

In the 1-D state, conservation of mass and conservation of energy place strong constraints on the dynamics (e.g. Ribe Reference Ribe1985). The total mass conservation equation, obtained by summing (4.1a) and (4.1b), can be readily integrated, subject to the impenetrable boundary condition ($u_z^s = u_z^l = 0$ at $z=0$) to give

(4.4)\begin{equation} \overline{\rho u}_z \equiv \lambda_\rho \phi^s \rho^s u_z^s + \phi^l \rho^l u_z^l = 0 ,\end{equation}

and thus the mean vertical velocity of the mixture $\bar {u}_z = 0$. Since $\phi ^\varepsilon, \rho ^\varepsilon$ are positive definite, the velocities of solid and liquid phases have to be of opposite sign. Physically, we expect solutions where dense (solid) phase falls ($u_z^s<0$) and light (liquid) phase rises ($u_z^l>0$). Further, integrating the solid mass conservation equation (4.1a) subject to $u_z^s = 0$ on $z = 0$ gives

(4.5)\begin{equation} \int _0^1 \varGamma^s \, \mathrm{d} z = \lambda_\rho ( \phi^s \rho^s u_z^s )|_{z=1}.\end{equation}

This relation states that the total melting (freezing) within the layer must balance with the solid mass flux entering (leaving) the layer. Note that the solid flux on the top boundary is not imposed through the boundary conditions (4.2); instead, it is determined as part of the solution. Since we expect the solid phase to always be falling, $u_z^s<0$, the integral of $\varGamma ^s$ is negative and thus on average heat is being absorbed by melting. Conservation of mass (4.4) also tells us that the solid mass flux in is balanced by the liquid mass flux out.

The total heat balance in the 1-D system, which follows from integrating the temperature equation (4.1e), may be written as

(4.6)\begin{equation} Q^C = Q^L + Q^P + Q^F + Q^V + Q^\zeta.\end{equation}

The term on the left-hand side represents the net conductive heat flux across the slurry layer, i.e. difference between heat flux out at the top $q_{z=1}^C$ and heat flux in at the bottom $q_{z=0}^C$:

(4.7)\begin{equation} Q^C = \int _0^1 - \frac{\mathrm{d} }{\mathrm{d} z} \left( \bar{k} \frac{\mathrm{d} T}{\mathrm{d} z} \right) \, \mathrm{d} z =-\left.\left( \bar{k} \frac{\mathrm{d} T}{\mathrm{d} z} \right)\right|_{z=1} + \bar{k} \theta = q_{z=1}^C - q_{z=0}^C. \end{equation}

Terms on the right-hand side represent contributions due to latent heat release/absorption $Q^L$, the heat-pipe effect $Q^P$, interphase friction $Q^F$, viscous heating $Q^V$ and heating due to compaction $Q^\zeta$:

(4.8)\begin{align} \left.\begin{array}{c@{}} \begin{aligned} Q^L & = \dfrac{\sqrt{R\, Pr}}{St} \int _0^1 \varGamma^s \, \mathrm{d} z,\\ Q^P & =-\dfrac{\lambda_\rho \sqrt{R\, Pr}}{St} \int _0^1 \dfrac{j_z }{ T} \dfrac{\mathrm{d} T}{\mathrm{d} z} \, \mathrm{d} z ,\\ Q^F & = D K \sqrt{R\, Pr}\int _0^1 j_z \Delta u_z \, \mathrm{d} z, \end{aligned} \quad \begin{aligned} Q^V & =\dfrac{4}{3} D Pr\int _0^1 \lambda_\mu \phi^s \left( \dfrac{\mathrm{d} u_z^s}{\mathrm{d} z} \right)^2 + \phi^l \left( \dfrac{\mathrm{d} u_z^l}{\mathrm{d} z} \right)^2 \, \mathrm{d} z,\\ Q^\zeta & = \lambda_\mu D Pr \int _0^1 \phi^s \mathcal{C}(\phi^l) ( \zeta^s )^2 \, \mathrm{d} z. \end{aligned} \end{array}\right\} \end{align}

(Recall that the condition (4.4) leads to elimination of terms involving $\bar {u}_z$ from the temperature equation (4.1e); thus the advective terms do not contribute to the heat balance.) We can deduce the sign of each contribution to the heat flux. The total melting rate balances the solid flux from above (4.5), and so $Q^L < 0$; heat is absorbed through melting. Since $u^s < 0$ and $u^l > 0$, the flux $j_z$ is always negative $j_z < 0$. The temperature gradient is also expected to be negative $\mathrm {d}{T}/\mathrm {d}{z}<0$ as temperature decreases monotonically with height. Thus, the product $j_z (\mathrm {d}{T} / \mathrm {d}{z}) > 0$, and since other quantities in the integrand are positive definite, it follows that the total contribution on account of the heat-pipe effect is negative $Q^P < 0$. Clearly, integrands corresponding to frictional, viscous and compaction heating are positive definite and thus provide positive contributions to the net heat flux $Q^F > 0, Q^V > 0\ {\rm and}\ Q^\zeta > 0$.

In a case where the mechanical heating terms ($Q^P, Q^F, Q^V, Q^\zeta$) are negligible, the main energy balance is between conductive heat flux and latent heat $Q^C = Q^L$; then, the total solid melting rate (integral of $\varGamma ^s$) depends only on the difference between heat flux in and out of the system (recall (4.7)). Note that the boundary conditions (4.2) only specify the heat flux at the bottom boundary; the heat flux out at the top boundary is freely determined by the solution. Thus, we do not know a priori what the difference is between heat flux in and out of the system, and consequently the total solid melting rate.

The pressure gradient is linked to the temperature gradient through the liquidus relation. From the imposed temperature and temperature gradient at the bottom boundary, the pressure gradient at $z=0$ is

(4.9)\begin{equation} \left.\frac{\mathrm{d} P}{\mathrm{d} z}\right|_{z=0} = \frac{-\theta}{( 1-\lambda_\rho^{-1} ) St D}.\end{equation}

We can compare this with the magnitude of the hydrostatic pressure gradient $\mathrm {d} P_H /\mathrm {d} z = - \bar {\rho }$. Neglecting density variations of individual phases (i.e. $\mathrm {d}\rho ^s=\mathrm {d}\rho ^l=0$), the mixture density is $1 \leq \bar {\rho } \leq \lambda _\rho$ (in dimensionless units). Therefore, the minimal pressure gradient required to support the weight of the matter is $-1$, and a crude bound on the pressure gradient imposed at $z=0$ is ${\mathrm {d} P}/{\mathrm {d} z}|_{z=0} < -1$. This, together with (4.9), provides a constraint for the required heat flux at $z=0$:

(4.10)\begin{equation} \theta > ( 1-\lambda_\rho^{-1} ) St D. \end{equation}

Clearly the parameters $\lambda _\rho, St$ and $D$ can significantly affect the nature of the basic state and play an important role in regulating the lower bound for the heat flux required to support the slurry layer. In the context of the slurry F-layer where $St = 4, D = 0.2$ and $\lambda _\rho = 1.02$, this bound requires $\theta > 0.016$, which is certainly compatible with geophysical estimates of $\theta$ (table 1).

Finally, note that we can use (4.4) to express $\phi ^s$ in terms of velocities and densities:

(4.11)\begin{equation} \phi^s = \frac{\rho^l u_z^l}{\rho^l u_z^l-\lambda_\rho \rho^s u_z^s } \quad \text{on } 0 < z \leq 1, \end{equation}

which provides a way to interpret volume fraction in terms of balance of momenta of the two phases. The above expression, however, could not be used to determine $\phi ^s$ on the impenetrable boundary ($z=0$) where both solid and liquid velocities are zero. Instead, the expression for $\phi ^s$ on the bottom boundary results from evaluation of the total mass balance (sum of (4.1a) and (4.1b)) at $z=0$, which yields $\phi ^s = \rho ^l ( {\mathrm {d} u_z^l}/{\mathrm {d} z} )/( \rho ^l ( {\mathrm {d} u_z^l}/{\mathrm {d} z} )-\lambda _\rho \rho ^s ( {\mathrm {d} u_z^s}/{\mathrm {d} z} ) )$.

4.3. Dominant balance and limiting behaviour

The two-phase model (3.7) is quite general in that it can be applied to a variety of scenarios involving fully compressible, inertial, two-phase flows. In particular, we have made no assumptions about the balance of forces and the nature of the flow. Here we derive a reduced 1-D model that will be helpful when interpreting results obtained by solving the more general 1-D equations. We initially consider the case without compaction viscosity $\mathcal {C}=0$, and come back to consider its effect at the end of the section.

We use two approximations that are commonly employed in studies of two-phase solid–liquid systems. One of them is the Boussinesq approximation which ignores density differences except where they appear in terms multiplied by the acceleration due to gravity $\boldsymbol { g }$. According to the Preliminary Reference Earth Model (PREM) (Dziewonski & Anderson Reference Dziewonski and Anderson1981), variations in core density are of the order of 0.1 % within hundreds of kilometres of the ICB. Thus, the Boussinesq approximation may be applied, whereby (in dimensionless units) $\rho ^s = \rho ^l = 1$ everywhere except in the buoyancy term (see § C.3). We also assume that both phases behave as very viscous fluids whose acceleration and inertia are negligible (see § C.4); thus terms ${\mathrm {D}^\varepsilon \boldsymbol {u}^\varepsilon }/{\mathrm {D} t}$ and $\varGamma ^\varepsilon \Delta \boldsymbol {u}$ are negligible, and the equations of solid and liquid momenta ((4.1c) and (4.1d)) simplify to

(4.12)$$\begin{gather} 0 =- \phi^s \frac{\mathrm{d} P}{\mathrm{d} z} - \lambda_\rho \phi^s \rho^s - Kj_z + \frac{4}{3} \lambda_\mu \sqrt{ \frac{Pr}{R}} \frac{\mathrm{d} }{\mathrm{d} z} \left( \phi^s \frac{\mathrm{d} u_z^s}{\mathrm{d} z} \right), \end{gather}$$

(4.13)$$\begin{gather}0 =- \phi^l \frac{\mathrm{d} P}{\mathrm{d} z} - \phi^l \rho^l + Kj_z + \frac{4}{3} \sqrt{ \frac{Pr}{R}} \frac{\mathrm{d} }{\mathrm{d} z} \left( \phi^l \frac{\mathrm{d} u_z^l}{\mathrm{d} z} \right). \end{gather}$$

On eliminating the pressure gradient term between (4.12) and (4.13), we obtain an expression for the solid flux:

(4.14)\begin{equation} \frac{Kj_z}{\phi^s \phi^l} =- ( \lambda_\rho \rho^s - \rho^l ) + \frac{4}{3} \sqrt{ \frac{Pr}{R}} \left\lbrace \frac{\lambda_\mu}{\phi^s} \frac{\mathrm{d} }{\mathrm{d} z} \left( \phi^s \frac{\mathrm{d} u_z^s}{\mathrm{d} z} \right) - \frac{1}{\phi^l} \frac{\mathrm{d} }{\mathrm{d} z} \left( \phi^l \frac{\mathrm{d} u_z^l}{\mathrm{d} z} \right) \right\rbrace .\end{equation}

The first term on the right-hand side of (4.14) represents the relative buoyancy between the two phases, whereas the second term represents the resistive effects of viscous stresses. We expect the dimensionless solid flux $j_z = \phi ^s \rho ^s \phi ^l \rho ^l \Delta u_z/\bar {\rho }$ to be maximised in the absence of viscous effects. This gives us an upper bound order of magnitude estimate for the relative velocity: ${\max | \Delta u_z | \sim ( \lambda _\rho -1 )/K}$ (using $\rho ^s = \rho ^l = \bar {\rho } = 1$). Moreover, since $u_z^s < 0$ and $u_z^l > 0, | \Delta u_z |= ( | u_z^s | + | u_z^l | )$, and thus the magnitude of the velocity of each phase is also of the order $\max | u_z^s | \sim \max | u_z^l | \sim ( \lambda _\rho -1 )/K$. In dimensional units, we have $\max | \Delta u_z | \approx (\lambda _\rho - 1) g /(\lambda _\rho c_D )$. Using values from table 1 yields an estimate of maximal relative velocities of the order of $5\times 10^{-5}$–$10^{-2}\ \textrm {m} \ \textrm {s}^{-1}$, which is much faster than the rate at which the inner core grows (a few millimetres per year, or $10^{-11} \textrm {m} \ \textrm {s}^{-1}$) and comparable to the estimated outer core flow speeds of $O(10^{-4})\ \textrm {m} \ \textrm {s}^{-1}$.

Since the viscosity of the solid phase is much greater than the viscosity of the liquid phase ($\lambda _\mu \gg 1$) we may neglect the viscous stress of the liquid phase in (4.14). Furthermore, the mass conservation constraint (4.4) allows us to write $u_z^s = -\phi ^l \rho ^l u_z^l/(\lambda _\rho \phi ^s \rho ^s)$, and $\Delta u_z = -\bar {\rho } u_z^l/(\lambda _\rho \phi ^s \rho ^s) = \bar {\rho } u_z^s/(\phi ^l \rho ^l)$. We can thus recast (4.14) as expressions for the phase velocities $u_z^l$ and $u_z^s$:

(4.15)$$\begin{gather} u_z^l = \frac{\lambda_\rho \phi^s}{ K \rho^l}( \lambda_\rho \rho^s - \rho^l ) - \frac{4}{3} \frac{\lambda_\mu}{K} \sqrt{ \frac{Pr}{R}} \frac{\lambda_\rho }{ \rho^l } \frac{\mathrm{d} }{\mathrm{d} z} \left( \phi^s \frac{\mathrm{d} u_z^s}{\mathrm{d} z} \right), \end{gather}$$

(4.16)$$\begin{gather}u_z^s =- \frac{ \phi^l}{ K \rho^s} ( \lambda_\rho \rho^s - \rho^l ) + \frac{4}{3} \frac{\lambda_\mu}{K} \sqrt{ \frac{Pr}{R}} \frac{1}{ \rho^s } \frac{\phi^l}{\phi^s} \frac{\mathrm{d} }{\mathrm{d} z} \left( \phi^s \frac{\mathrm{d} u_z^s}{\mathrm{d} z} \right). \end{gather}$$

It is useful to note that the term under the $z$-derivative in (4.15) may be rewritten with use of (4.4) as

(4.17)\begin{equation} \phi^s \frac{\mathrm{d} u_z^s}{\mathrm{d} z} = \frac{\phi^l}{\lambda_\rho} \left( ( \lambda_\rho u_z^s - u_z^l )\frac{1}{\phi^l}\frac{\mathrm{d} \phi^l}{\mathrm{d} z} - \frac{\mathrm{d} u^l}{\mathrm{d} z} \right). \end{equation}

This is useful because it reveals that the viscous stress term is proportional to $\phi ^l$ (rather than to $\phi ^s$ as could be interpreted at first glance from the left-hand side of (4.17)).

We can now estimate how $u^s, u^l$ and $\phi ^s$ ($\phi ^l$) behave when parameters $R, \lambda _\mu$ and $K$ are very large. To do so, we assume that the two terms which set the velocities of the two phases (relative buoyancy and viscous resistance) in both (4.15) and (4.16) have to remain in balance as parameters $R, \lambda _\mu$ and $K$ are increased to large, but finite, values. We also assume that the $z$-derivatives are of order unity throughout the bulk of the domain, which is consistent with the scaling estimates we infer.

When $R \gg 1$, the coefficient of the second term in (4.15) becomes very small. For the two terms in (4.15) to remain in balance would require the buoyancy term to similarly decrease, and thus $\phi ^s \sim R^{-1/2}$. It follows then that $u^l$ should also scale with $R^{-1/2}$. Since by mass conservation (4.4) $u^s_z \sim {( 1 - \phi ^s ) u_z^l}/{\phi ^s}$, we expect $u_z^s$ to tend to become $R$-independent for $R \gg 1$. This is also consistent with the expression (4.16) (together with (4.17)), where the $R$ dependence cancels out from the viscous term if $\phi ^s \sim R^{-1/2}$. If $u^l \sim R^{-1/2}$, then the liquid-phase velocity decreases towards zero for as $R\to \infty$ and thus all of the solid flux $j_z$ has to be carried by the solid-phase velocity (recall (4.14)) and thus we expect that $u_z^s$ approximately tends to the ‘terminal velocity’: $-(\lambda _\rho -1)/K$. Since $j_z \sim \phi ^s \phi ^l \Delta u_z$, we expect the magnitude of the solid flux itself to decrease for large $R$; in the limit $R\to \infty, \phi ^s \sim R^{-1/2}$ ($\phi ^l \to 1$) and $\Delta u_z$ tends to a constant, and thus we expect that $j_z \sim R^{-1/2}$.

When $\lambda _\mu \gg 1$, the coefficient of the second term in (4.16) becomes very large. The two terms in (4.16) can stay in balance if $\phi ^l \sim \lambda _\mu ^{-1}$ (recall again (4.17), whereby the first term in (4.16) goes like $\phi ^l$, while the second goes like $\lambda _\mu {\phi ^l}^2$). It follows then that $u_z^s \sim \lambda _\mu ^{-1}$. Moreover, since $u_z^l \sim \phi ^s u_z^s /\phi ^l$, it is expected that $u_z^l$ becomes independent of $\lambda _\mu$ for $\lambda _\mu \gg 1$. As $u_z^s$ tends to zero for large $\lambda _\mu$, we expect $u_z^l$ to tend towards the ‘terminal velocity’ $\sim (\lambda _\rho -1)/K$. In the limit $\lambda _\mu \to \infty$, the solid flux decreases $j_z \sim \phi ^s\phi ^l \Delta u_z \sim \lambda _\mu ^{-1}$.

Using the same approximations and reasoning as above it is clear that the behaviour at small $R$ (large $R^{-1/2}$) is analogous to the behaviour at large $\lambda _\mu$, and so we have $u_z^s \sim \phi ^l \sim j_z \sim R^{1/2}$ and $u_z^l \sim (\lambda _\rho -1)/K$ for small $R$. Note, however, that the converse is not true: the small $\lambda _\mu$ limit is not the same as large $R$ limit. Once $\lambda _\mu$ becomes small, or of order unity, the neglect of the viscous term corresponding to the liquid phase in (4.14) is unjustified.

Finally, between the expressions (4.14)–(4.16) it is clear that the velocities of both phases (as well as the solid flux) should scale with $K^{-1}$ for large $K$, and we expect $u_z^s\sim u_z^l \sim j_z \sim K^{-1}$ as $K \rightarrow \infty$. The scaling behaviours of $u^s, u^l$ and $\phi ^s$ ($\phi ^l$) with respect to parameters $R, \lambda _\mu$ and $K$ are summarised in table 3.

Table 3. Summary of scalings.

We can also estimate the dominant balance in the heat equation. Since the phase velocities are at most of the order of $(\lambda _\rho -1)/K$ (small), we expect terms involving velocities in the the heat equation (4.1e) to be negligible. Thus, the leading-order balance in the heat equation (4.1e) is between latent heat and the conductive heat flux:

(4.18)\begin{equation} - \frac{\varGamma^s}{St} \approx \frac{1}{\sqrt{R\, Pr}} \frac{\mathrm{d} }{\mathrm{d} z} \left( \bar{k} \frac{\mathrm{d} T}{\mathrm{d} z} \right). \end{equation}

Finally, we discuss the effect of compaction. When density variations are small we have approximately $\zeta ^s = \mathrm {d} u_z^s/\mathrm {d} z$. Since $u_z^s < 0$, with $u_z^s = 0$ on the bottom boundary, we expect that $\mathrm {d} u_z^s/\mathrm {d} z < 0$ at least near the base of the domain. At the top boundary we have $\mathrm {d} u_z^s/\mathrm {d} z = 0$, so it is likely that $\mathrm {d} u_z^s/\mathrm {d} z < 0$ throughout the layer, and so $\zeta ^s < 0$, i.e. the solid phase is compacting and $\Delta p > 0$. All compaction terms appear multiplied by $\mathcal {C}(\phi ^l)$. Thus, effects of compaction are maximised in parts of the domain that are predominantly solid. Since we are mainly interested in slurry-type solutions where $\phi ^s$ is small in the bulk of the domain we expect the effects of compaction to be greatest near the lower boundary.

Compaction terms appear in the solid momentum equation, the liquidus relation and the heat equation. We can estimate the importance of compaction terms relative to other terms in these equations. The compaction term in the temperature equation is similar in nature to the viscous heating term, and similarly is of the order of $((\lambda _\rho -1)/K)^2$ (small squared). Of course, the $\mathcal {C}(\phi ^l)$ factor can exaggerate the magnitude of compaction heating in the portion of the domain with $\phi ^l\ll 1$. Nevertheless, this term is ultimately of the same order of magnitude as the other heating terms (frictional, viscous, heat-pipe) that are negligible to leading order. In the liquidus relation (4.1f) compaction terms enter paired with the density gradient. Not only is $\zeta ^s \sim (\lambda _\rho -1)/K$ small, but it is also multiplied by the small density gradient. Thus the effect of compaction is negligible in comparison with the other terms in that equation. We therefore expect that the primary influence of compaction terms is in the solid momentum equation.