A two-phase pure slurry model for planetary cores: one-dimensional solutions and implications for Earth's F-layer

Abstract

We develop and analyse a continuum model of two-phase slurry dynamics for planetary cores. Mixed solid–liquid slurry regions may be ubiquitous in the upper cores of small terrestrial bodies and have also been invoked to explain anomalous seismic structure in the F-layer at the base of Earth's liquid iron core. These layers are expected to influence the dynamics and evolution of planetary cores, including their capacity to generate global magnetic fields; however, to date, models of two-phase regions in planetary cores have largely ignored the complex fluid dynamics that arises from interactions between phases. As an initial application of our model, and to focus on fundamental fluid dynamical processes, we consider a non-rotating and non-magnetic slurry comprised of a single chemical component with a temperature that is tied to the liquidus. We study one-dimensional solutions in a configuration set up to mimic Earth's F-layer, varying gravitational strength $R$, the solid/liquid viscosity ratio $\lambda _{\mu }$ and the interaction parameter $K$, which measures friction between the phases. We develop scalings describing behaviour in the limit $R \gg 1$ and $\lambda _{\mu } \gg 1$, which are in excellent agreement with our numerical results. Application to Earth's core, where $R \sim 10^{28}$ and $\lambda _{\mu } \sim 10^{22}$, suggests that a pure iron slurry F-layer would contain a mean solid fraction of at most 5 %.

1. Introduction

Multiphase flows involving phase changes are ubiquitous in nature. Everyday examples include formation of clouds, rain droplets and snow in the atmosphere (Pruppacher 1986). In the Earth's mantle, two-phase flows occur in the form of melting of rock and subsequent transport of magma to the surface (e.g. McKenzie 1984; Spiegelman 1993; Bercovici, Ricard & Schubert 2001). Two-phase flows also arise in deep planetary interiors. Planetary differentiation driven by heating from radioactive isotope decay and planetary accretion leads to partial melting, and ultimately to the segregation of a molten iron-rich core from the solid silicate mantle (Rubie, Nimmo & Melosh 2015). Subsequent cooling leads to crystallisation of the fluid core, which drives compositional convection that provides an efficient power source for maintaining planetary magnetic fields by dynamo action (Breuer, Rueckriemen & Spohn 2015). In Earth's core, freezing of dense iron-rich material proceeds from the inside out, leading to segregation of solid and liquid phases with power provided by the continual enrichment of the liquid in light impurities (Nimmo 2015). By contrast, small planetary bodies such as Ganymede, the moon, and asteroids could freeze dense material from the top, with iron ‘snowing’ down from the core–mantle boundary (Breuer et al. 2015). In this case, power for dynamo action is thought to be provided both by the falling snow and by the gravitational energy released when the snow remelts at greater depths (Rückriemen, Breuer & Spohn 2015; Davies & Pommier 2018).

This work is motivated by another ‘snow zone’ that is hypothesised to exist at the base of Earth's liquid core. The traditional picture of the geodynamo is that turbulent convection in the outer core produces a state in which the interior is nearly homogeneous and adiabatic. However, seismic studies indicate the existence of a region at the base of the outer core where density is stably stratified (see Souriau & Poupinet (1991), and also the overview in Gubbins, Masters & Nimmo (2008)). The depth of this region, called the F-layer, is estimated to be between 150 and 400 km. The presence of such a region at the base of the outer core has potentially far-reaching consequences since the F-layer mediates the power input to the bulk core where the magnetic field is generated. It is therefore important to understand the extent of the influence that the F-layer has on the motions in the turbulent bulk of the core, and how that in turn influences the operation of the geodynamo.

Several possibilities have been proposed to explain the anomalous stratification of the F-layer. The simplest explanation, a thermal boundary layer, is not viable because the thermal gradient at the base of the liquid core is strongly destabilising (Gubbins et al. 2008). One possibility is that the layer results from convective translation of the inner core, a deformationless mode of motion that can arise if the inner core is convectively unstable. Translation results in the inner core freezing in the western hemisphere and melting in the eastern hemisphere, which produces an iron-rich dense layer that sits above the inner core boundary (ICB) (Alboussiere, Deguen & Melzani 2010). However, recent upward revisions of the thermal conductivity of iron (Pozzo et al. 2014) make thermally driven convective translation unlikely. Compositional (Gubbins, Alfe & Davies 2013) or double-diffusive (Deguen, Alboussière & Labrosse 2018) convection in the inner core can arise, but the conditions for these instabilities do not appear to be satisfied at the present day. It is also unclear whether convective translation can produce the observed global F-layer thickness. Gubbins et al. (2008) proposed a hydrostatic thermochemical model of iron alloyed with light element at the liquidus temperature that succeeded in producing a thick stably stratified F-layer. However, the stabilising compositional gradient that produces the overall density stratification is imposed through the boundary conditions and so its origin, together with the mechanism by which light material passes through the stratified layer, is not explained by the model.

An enticing hypothesis inferred by analogy with metallurgical studies of solidifying mixtures is that the region of the outer core close to the ICB is in a state of constitutional supercooling, and thus matter in the F-layer could be composed of a mixture of liquid and solid phases (Loper & Roberts 1978, 1981). Within the two-phase hypothesis, two distinct regimes were proposed: a slurry zone with a low solid fraction (e.g. Loper & Roberts 1978) and a mush zone with a high solid fraction (e.g. Fearn, Loper & Roberts 1981). In a slurry regime, solid phase exists in the form of fine grains suspended in the parent liquid phase. In a mush regime, solidification is assumed to form a dendritic solid skeleton with liquid flowing along pore spaces.

To explain the F-layer as a two-phase region, its properties must be consistent with seismic observations. The most robust observation is a flattening of the P-wave velocity profile with depth, which is usually interpreted to represent density stratification (Gubbins et al. 2008). This behaviour has been reproduced by simple thermodynamic slurry models (Wong, Davies & Jones 2021), but has not so far been reported for the mush regime. There is also no clear seismic evidence of shear-wave velocities, reflections or impedance contrasts above the ICB (Souriau 2007), which is consistent with the existence of a slurry. A mush can in principle support shear waves and top-side reflections owing to its high solid fraction but the seismic properties of a mushy layer above the ICB are unknown; so while these observations do not rule out a mushy F-layer they require that it is somehow hidden from current seismic probes, which is a solution we do not favour. The mush hypothesis does have the advantage that it provides natural nucleation sites, whereas nucleation in a pure slurry may require very substantial undercooling (e.g. Huguet et al. 2018). However, the nucleation difficulty, found in theoretical calculations for a pure substance, may not apply in the presence of impurities which may well be present in the core. Our view is that core nucleation is very imperfectly understood, and therefore it is not appropriate to rule out the slurry model on that basis. The solid fraction in a two-phase F-layer is unknown, but estimates based on simple thermodynamic arguments suggest values of the order of 1 % (Gubbins et al. 2008), while estimations based on mineral physics data give ${\sim }15\,\%$ (Zhang et al. 2019). These values are smaller than the rheological transition above which the two-phase region is expected to form an interconnected solid network (Solomatov 2007). Finally, a mushy zone at the base of Earth's core is expected to have a maximum thickness of a few kilometres, much thinner than the F-layer, since a thick mush would collapse under its own weight (Deguen, Alboussière & Brito 2007). Based on these considerations we focus on a slurry F-layer in the rest of this paper. In this scenario, heavy solids fall under gravity, while the light liquid rises buoyantly upwards. A net inward transport of dense solid and a net outward transport of light material can result in an overall stable stratification.

Previous models of the slurry F-layer (and indeed ‘iron snow’ models for planetary cores in general) have focused on the thermodynamics of the problem rather than the fluid dynamics. Wong, Davies & Jones (2018) and Wong et al. (2021) applied the two-phase ‘diffusive mixture’ model of a slurry developed by Loper & Roberts (1978) to study the F-layer and found solutions with density and P-wave velocity variations that are in good agreement with those inferred from seismic observations. The model of Wong et al. (2018) is purely thermodynamic in that it assumes that the F-layer is in hydrostatic equilibrium and on the liquidus, where the liquid phase is at rest and the sedimentation of solid phase is approximated by Stokes flow along the pressure gradient. The diffusive mixture formulation models the evolution of the mean mixture, with a single momentum equation for the barycentric flow velocity, so flow velocities of individual constituent phases cannot be known. Here we follow an alternative ‘two-fluid’ formulation which models each individual phase as a separate fluid continuum (Roberts & Loper 1987). In this more general model, which can be simplified to the diffusive model, the velocities of each phase as well as the solid sedimentation flux are determined self-consistently. Our study therefore represents the first steps towards a fully fluid dynamical model of the two-phase slurry F-layer.

Throughout this work we consider a pure iron slurry, i.e. one composed of a single chemical constituent. This is a simplification of the real core composition, which is expected to comprise one or more light elements. However, we will show that the pure slurry is a rich and complex fluid dynamical system that is worthy of study in its own right. Moreover, by elucidating the behaviour of this system across a wide range of control parameters we obtain a crucial reference point for understanding the additional complexities due to adding lighter elements, which will be considered in future work. Finally, the pure system is of geophysical interest because it allows us to test previous predictions based solely on thermodynamic arguments (Gubbins et al. 2008) that any slurry that could form in the F-layer would contain only a small fraction of solid phase.

Modelling multiphase dynamics of solid–liquid slurry poses a significant theoretical challenge. Broadly speaking, there are two strategies for treating the solid particles: as rigid entities or as a continuum field. A basic example of the rigid approach is modelling of individual sinking solids and their experience of drag and their effect on the ambient chemical gradients (e.g. Inman et al. 2020; Magnaudet & Mercier 2020). This strategy can be extended to models involving many tracers, representing solid crystals, in a convective flow (e.g. Suckale et al. 2012; Patočka, Calzavarini & Tosi 2020). However, computational expense limits the number of tracers that can be modelled, while the addition, removal or change of sizes of particles during simulation as well as the associated latent heat exchanges are challenging to represent numerically. With this approach, numerical modelling at the scale of the F-layer is currently inaccessible.

When the solid phase is treated as a continuum, macroscopic behaviour resulting from microscopic phase interactions is modelled by interpenetrating and interacting continuum fields. The exact locations of the individual phases, and hence the interfaces between phases and the microscopic properties of melting, are unknown; at the system scale, phases are represented by the volume fraction they occupy in a control volume. This approach therefore lacks a first-principles link to local-scale variables and processes (e.g. formation and growth of solid crystals or topology and evolution of phase interfaces) and thus requires that local mechanical and thermodynamic phase interactions are parametrised through constitutive relations based on macroscopic variables. The deliberate focus on the system scale constitutes both the key limitation but also the chief strength of continuum models. By abandoning a rigorous representation of local-scale phase interactions, the continuum framework allows one to construct models of system-scale behaviour of tractable mathematical complexity. Such models have a long history (e.g. Truesdell & Toupin 1960) and have had significant impact in geophysics, as highlighted by the two-phase model of Drew (1971) and Drew & Segel (1971) and its subsequent application to the problem of two-phase melt transport by McKenzie (1984) and others (e.g. Fowler 1985; Ribe 1985; Bercovici et al. 2001). In recent years, models of this kind have found many areas of application, such as magma oceans (Boukaré & Ricard 2017), planetary mantle evolution (Zhang, Bercovici & Jordan 2021), ice–water transport on Europa (Kalousova et al. 2014) and magnetic induction in mushy zones (Bercovici & Mulyukova 2021).

Within the two-phase continuum theory, a number of new phenomena arise that are not present in a single-phase flow. In a pure slurry, in addition to the effects of thermal buoyancy, convection of the liquid can be driven by buoyancy of solid grains. Furthermore, phase change and associated absorption/release of latent heat introduce an additional avenue of heat transport through the system. This gives rise to the ‘heat-pipe’ effect, whereby heat is transported by fluid parcels that solidify at their destination, releasing latent heat.

In this paper we study two-phase dynamics of a pure iron slurry in an effort to shed light on the dynamics in the F-layer. We extend the formulation of Roberts & Loper (1987) and derive and solve a set of equations that describe the motion of a reactive two-phase system. As described above, the bulk of the outer core (above the F-layer) is vigorously convecting, which mixes this region to almost uniform composition and an adiabatic temperature profile. On the other hand, the F-layer is stable to convection so we expect that any motion in the slurry will be more temperate, consisting chiefly of settling of solid phase under gravity. We further expect that radial variations will dominate the signals obtained by seismology because the inferred density anomalies are many orders of magnitude larger than the thermochemical anomalies that drive core convection (Stevenson 1987). Hence, in this paper, we restrict our attention to one-dimensional (1-D) solutions representing motion along the vertical direction, which serve as a natural starting point from which to build more detailed dynamical models. We are interested in gaining a broad understanding of the nature of solutions to these two-phase equations up to solid fractions of $\phi ^s \sim 0.5$, beyond which a mushy zone is expected to form. In particular, we seek to understand the factors that determine the distribution of solid phase within the layer, and to provide independent constraints on the solid fraction that can exist in a pure slurry at the physical conditions of Earth's F-layer.

Previous applications of two-phase modelling to planetary mantle dynamics have mainly focused on situations where the solid phase is dominant (e.g. McKenzie 1984; Bercovici et al. 2001; Šrámek, Ricard & Bercovici 2007; Boukaré & Ricard 2017; Zhang et al. 2021). These studies assume that the liquid volume fraction stays low and the solid phase forms an interconnected matrix. In such a case inertia of both phases is negligible and momentum balance for the liquid phase is reduced to Darcy's law, whereby the liquid percolation through the solid is determined by the overpressure and the permeability of the solid matrix. By contrast, in the slurry limit the solid phase is expected to be sparse. In this case it is unclear a priori which terms in the equations can be neglected. In this study we take a general approach that does not assume Darcy flow or negligible kinetic energy. We draw heavily from the two-phase equations for a pure slurry developed by Roberts & Loper (1987), which, to our knowledge, have not previously been solved. We attempt to strike a balance between simplicity and realism, noting that many properties of deep Earth materials are poorly known owing to the immense pressures and temperatures that arise in planetary cores.

The outline of the paper is as follows. In § 2 we describe the governing equations for the two-phase slurry. In § 3 we formulate the mathematical problem of the two-phase slurry in the context of the Earth's F-layer. In §§ 4 and 5 we analyse and solve the equations for the case of 1-D time-independent solution. Concluding discussion is presented in § 6.

2. Equations governing two-phase flow

Derivation of conservation laws in two-phase theory involves averaging of the microscopic scale equations to obtain equations governing the macroscopic scale (see e.g. Drew 1971, 1983). The averaging formalism for deriving continuum equations of macroscopic two-phase flow is well established, and does not need to be reproduced here. In this section we describe the continuum equations governing the evolution of a two-phase slurry. Our approach follows closely the work of Drew (1971, 1983) and Roberts & Loper (1987) and also uses results from Bercovici et al. (2001) and Šrámek et al. (2007). The complete set of equations formulated in this section may be found in Appendix C. A table of symbols involved in the development of the governing equations in this section is provided in Appendix A (table 4). Readers more interested in the analysis of the mathematical problem may pass directly to § 3.

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. 2001) 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 1983; Roberts & Loper 1987; Bercovici et al. 2001).

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 (1987) 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 1978).

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 1983). 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 (2010), 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 (1987) and Bercovici et al. (2001) 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 2003). 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 (1987), 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 1987).

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 (2007) (see also Rudge, Bercovici & Spiegelman 2011), 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 2007) (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 1987), 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 1987; Ricard 2007; De Groot & Mazur 2013) 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 1988), and acts to reduce the velocity differences between phases. In the tight-coupling limit $\varLambda _D \to \infty$ (Roberts & Loper 1987), 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. 2001; Michaut et al. 2013). 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 2010; Degruyter et al. 2012; Michaut et al. 2013). 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 2019). 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 1987). Using the maximum grain size of 3 cm as estimated by Walker, Davies & Wilson (2021), 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 1987).

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. 2022).

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 2003). 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 2007). Bercovici & Ricard (2003) 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 (1987) 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 1987). 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 (2007) and Loper & Roberts (1978) 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 2007). 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 2007), 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 1978):

(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. Application to the F-layer

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 1984).

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. (2021). 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 2020) to $5\times 10^{-10}$ (Šrámek et al. 2007). 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. 2008). 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. One-dimensional time-independent solution

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 1985). 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 1981), 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.

5. Results

In this section we present the solutions to the full system of (4.1), and describe the effects of varying $R, \lambda _\mu$ and $K$ in turn. For now, we set $\mathcal {C}=0$ and return to investigate the effects of compaction in § 5.6.

5.1. Variation of $R$

Figure 2 shows profiles of the freezing rate $\varGamma ^s$, the solid volume fraction $\phi ^s$ and vertical velocities of solid and liquid phases $u_z^s, u_z^l$, for increasing values of $R$ (at fixed $\lambda _\mu =10^5$ and $K=10^3$). (Our discussion is focused on the variables plotted in figure 2; for completeness plots of density, temperature and pressure gradient are included in Appendix E, figure 13.) Increasing $R$ enhances the downward gravitational force and causes the velocity of falling solid phase $u_z^s$ to increase in magnitude (become more negative) and that of rising liquid phase $u_z^l$ to decrease. A natural consequence of larger (more negative) $u_z^s$ and smaller $u_z^l$ is a reduction of solid volume fraction $\phi ^s$ throughout the layer (recall (4.11)). Solid material falls through the layer and converges ($\phi ^s \rightarrow 1$) near the bottom boundary. As $R$ increases, the gradients of $\phi ^s$ and $u_z^s$ near the boundary become larger and this convergence of solid mass becomes magnified. As a result of mass conservation (4.4) we observe a surge in the liquid velocity $u^l_z$ near the bottom boundary, where $\phi ^l \ll 1$, and $u^l_z \sim {\phi ^s u_z^s}/{\phi ^l}$ must be very large (locally) in order to balance the mass flux. By mass conservation (4.1a) the convergence of solid mass is balanced by melting; thus increasing $R$ results in more melting near the bottom, on account of increased convergence of solid mass.

Figure 2. Profiles of freezing rate $\varGamma ^s$ (a), solid volume fraction $\phi ^s$ (b), solid-phase velocity $u_z^s$ (c) and liquid-phase velocity $u_z^l$ (d) for increasing values of $R$; other parameters are fixed at $K = 10^3, \lambda _\mu = 10^5$. Recall that negative values of $\varGamma ^s$ indicate melting.

It is useful to define an average of any scalar quantity $A$ as

(5.1)\begin{equation} \langle A \rangle = \int _0^1 A \, \mathrm{d} z. \end{equation}

Figure 3 shows the variation of average solid velocity, average liquid velocity and average solid volume fraction versus $R$. Average solid velocity increases with $R$, exhibiting the expected $R^{1/2}$ behaviour at low $R$. The liquid velocity decreases with $R$ and, at large $R$, follows the scaling $R^{-1/2}$ (figure 3b), which is consistent with the dominant balance prediction of § 4.3. As $R$ is increased, average solid volume fraction decreases in line with dominant balance scalings. For large $R$ we observe that $\langle \phi ^s \rangle \sim R^{-1/2}$ (as figure 3c), and for small $R$ we observe that $\langle \phi ^l \rangle = 1 - \langle \phi ^s \rangle \sim R^{1/2}$ (see Appendix E, figure 14a). Since $\langle u_z^l \rangle \sim \langle \phi ^s \rangle \sim R^{-1/2}$ for large $R$, we expect $\langle u_z^s \rangle$ to become independent of $R$ for $R \gg 1$ (§ 4.3) and there is evidence of this behaviour in figure 3(a). The scaling of $u_z^s \sim R^{1/2}$ together with $\phi ^l \sim R^{1/2}$ implies that $\langle u_z^l \rangle$ becomes independent of $R$ at low $R$ and again this is seen in figure 3(b). For completeness we include the plot of the variation of the average solid flux in Appendix E (figure 14b), which verifies the predicted limiting behaviours: $\langle \,j_z \rangle \sim R^{1/2}$ as $R \to 0$ and $\langle \,j_z \rangle \sim R^{-1/2}$ as $R \to \infty$. Overall we find that the reduced model developed in § 4.3 captures the gross behaviour of these solutions.

Figure 3. Variation of the average solid-phase velocity $K|\langle u_z^s \rangle | / (\lambda _\rho -1)$ (a), the average liquid-phase velocity $K\langle u_z^l \rangle / (\lambda _\rho -1)$ (b) and the average solid volume fraction $\langle \phi ^s \rangle$ (c) with respect to $R$. In (a,b) the phase velocities have been scaled by the estimated terminal velocity $(\lambda _\rho - 1)/K$.

5.2. Variation of $\lambda _\mu$

Figure 4 shows solutions for increasing value of the viscosity ratio $\lambda _\mu$. As $\lambda _\mu$ is increased the sharp vertical variation of $u_z^s$ and $\phi ^s$ near the bottom boundary becomes completely smoothed out. The concurrent decrease in gradients of $u_z^s$ and $\phi ^s$ means reduced convergence of the solid phase and thus reduced melting (recall (4.1a)). Liquid-phase velocity increases throughout the bulk of the layer. The sharp variation of $u_z^l$ near $z=0$ remains unaffected.

Figure 4. Profiles of freezing rate $\varGamma ^s$ (a), solid volume fraction $\phi ^s$ (b), solid-phase velocity $u_z^s$ (c) and liquid-phase velocity $u_z^l$ (d) for increasing values of the viscosity ratio $\lambda _\mu$; other parameters are fixed at $K = 10^3, R=10^8$.

Variation of average quantities is presented in figure 5. As $\lambda _\mu$ is increased, viscous forces act to reduce $\langle u_z^s \rangle$ (figure 5a) and the average liquid content $\langle \phi ^l \rangle = 1 - \langle \phi ^s \rangle$ decreases. For large $\lambda _\mu$ the change in $\langle u_z^s \rangle$ and $\langle \phi ^l \rangle$ is proportional to $\lambda _\mu ^{-1}$, as found in § 4.3. The average liquid-phase velocity $\langle u_z^l \rangle$ increases with $\lambda _\mu$, but becomes independent of $\lambda _\mu$ in the limit $\lambda _\mu \gg 1$. This is expected since $u_z^l \sim \phi ^s u_z^s/\phi ^l$, with $u_z^s \sim \phi ^l \sim \lambda _\mu ^{-1}$ for $\lambda _\mu \gg 1$, the $\lambda _\mu$ dependence cancels out. The overall behaviour is consistent with our expectation from the reduced model (recall table 3).

Figure 5. Variation of the average solid-phase velocity $K|\langle u_z^s \rangle | / (\lambda _\rho -1)$ (a), the average liquid-phase velocity $K\langle u_z^l \rangle / (\lambda _\rho -1)$ (b) and the average liquid volume fraction $\langle \phi ^l \rangle$ (c) with respect to the viscosity ratio $\lambda _\mu$. In (a,b) the phase velocities have been scaled by the estimated terminal velocity $(\lambda _\rho - 1)/K$.

5.3. Variation of $K$

The interphase friction force acts to reduce the relative velocity between the two phases. In fact, the limit $K\to \infty$ is usually described as the tight-coupling limit, where the two phases move with the same velocity ($\Delta \boldsymbol {u} = 0$). We expect the same trend to apply here. Since $u_z^s < 0$ and $u_z^l>0$ always, reduction in the relative velocity can only be accomplished by an increase in $u_z^s$ (slower falling of the solid phase) and a decrease in $u_z^l$ (slower rise of the liquid phase).

Figure 6 shows solution profiles for increasing values of $K$. As $K$ is increased the flow velocities of both solid and liquid phases are reduced, and their vertical gradients near the bottom boundary become smoothed out due to stifled flow. The resulting flow velocities at large $K$, with $u_z^s$ nearly uniform and $u_z^l$ gradually increasing with height, produce a peculiar distribution of solid volume fraction, whereby solid content in the interior of the layer is smaller than the solid content on either boundary. Figure 7 shows that as $K$ is increased, average velocities of both phases decrease. For sufficiently large $K$, the decrease in velocities follows scaling $u_z^s \sim u_z^l \sim K^{-1}$. Average solid fraction also increases with $K$, but becomes independent of $K$ in the limit $K \to \infty$ (this, again, is consistent with the mass conservation constraint; see (4.11)).

Figure 6. Profiles of freezing rate $\varGamma ^s$ (a), solid volume fraction $\phi ^s$ (b), solid-phase velocity $u_z^s$ (c) and liquid-phase velocity $u_z^l$ (d) for increasing values of the interphase friction parameter $K$; other parameters are fixed at $R=10^8, \lambda _\mu = 2 \times 10^5$.

Figure 7. Variation of the average solid-phase velocity $|\langle u_z^s\rangle |$ (a), the average liquid-phase velocity $|\langle u_z^l\rangle |$ (b) and the average solid volume fraction $\langle \phi ^s\rangle$ (c) with respect to the interphase friction parameter $K$.

5.4. Regime diagram

In the above sections we encountered two distinct types of limiting behaviour: a stagnant solid regime characterised by large solid volume fraction and small solid-phase velocity; and a stagnant liquid regime characterised by large liquid volume fraction and small liquid-phase velocity. The stagnant solid limit is accessed as $R\to 0$, or $\lambda _\mu \to \infty$, whereas the stagnant liquid limit is reached as $R \to \infty$. These limits are visualised in figure 8, which shows regime diagrams of the average solid volume fraction and average phase velocities (normalised by $(\lambda _\rho -1)/K$) in the $R$–$\lambda _\mu$ space for three values of the interaction parameter $K = (10^2, 10^3, 10^4)$.

Figure 8. From top to bottom: contours of the average solid volume fraction $\langle \phi ^s \rangle$ (a–c); (normalised) average velocity of the solid phase $K | \langle u_z^s \rangle |/(\lambda _\rho -1)$ (d–f); (normalised) average velocity of the liquid phase $K \langle u_z^l \rangle /(\lambda _\rho -1)$ (g–i). Values of $K$ increase from left to right: (a,d,g) $K = 10^2$; (b,e,h) $K = 10^3$; (c,f,i) $K = 10^4$. Empty regions in the contour plots correspond to parameter values for which we do not have solutions.

For all three variables ($\langle \phi ^s \rangle, \langle u^s \rangle, \langle u^l \rangle$), the largest change in the slopes of the contour plots as $K$ increases occurs in region where $K$ is not much smaller that either $R$ or $\lambda _\mu$. From the plot, it appears that the largest change in the morphology of contour lines occurs in the regions where $R \lesssim K^2$ and $\lambda _\mu \lesssim K^2$. On the other hand, when both $R$ and $\lambda _\mu$ are much larger than $K$, changing $K$ has a more moderate effect on the characteristics of the solution. Indeed, we observe that regions where $R \gtrsim K^2$ and $\lambda _\mu \gtrsim K^2$ are remarkably self-similar to one another in that the contours share similar morphology and inclination, and only appear shifted across $R$–$\lambda _\mu$ parameter space. This self-similarity may be better visualised by zooming in on the regions $R>K^2, \lambda _\mu >K^2$, which we provide in Appendix E (figure 15).

Contours of $\langle \phi ^s \rangle$ in figure 8(a–c) run diagonally across the $R$–$\lambda _\mu$ parameter space for large $R$ and $\lambda _\mu$ (see also figure 15). The transition between the ‘stagnant solid’ (low-$R$) and ‘stagnant liquid’ (large-$R$) behaviour takes place roughly when $R \approx \lambda _\mu$ (see for example 9a). This information, together with the inferred behaviour of $\langle \phi ^s \rangle$ at both limits, allows extrapolation to parameter values relevant to planetary cores.

Figure 9. (a) Variation of average solid fraction $\langle \phi ^s \rangle$ together with lines of fit based on the scaling of $\phi ^s$ as $R\to 0$ (dashed lines) and as $R\to \infty$ (dash-dotted lines). Dotted lines indicate where the two lines of fit cross. (b) Extrapolation of average solid fraction $\langle \phi ^s \rangle$ to large values of $R$ and $\lambda _\mu$. Red rectangle highlights the estimated range of $R$ and $\lambda _\mu$ for the F-layer (table 2).

Figure 9(b) shows an extrapolation of the average solid volume fraction to planetary values of $R$ and $\lambda _\mu$. The viscosity contrast $\lambda _\mu$ is the parameter with the largest range of uncertainty, with estimates ranging from $10^{13}$ to $10^{24}$ (recall tables 1 and 2). Any slurry that could form in the F-layer would contain a very small fraction of solid. Gubbins et al. (2008) estimated a solid fraction of ${\sim }1$ %, which is consistent with our results. The low $\langle \phi ^s \rangle$ puts a constraint on the allowable combinations of $\lambda _\mu$ and $R$. Assuming $\langle \phi ^s \rangle \lessapprox 0.01$ in the F-layer requires the viscosity contrast $\lambda _\mu$ to be significantly smaller than the gravitational strength $R$. Taking $R = 10^{27}$ would require $\lambda _\mu \lesssim 10^{22}$, which is consistent with mineralogical determinations.

5.5. Heat flux

Figure 10 shows contributions of each term in the heat balance (4.6) as a function of $R, \lambda _\mu$ and $K$. Clearly, the contribution due to latent heat dominates at all parameter values considered. Thus, to leading order, the heat balance (4.6) is

(5.2)\begin{equation} q_{z=1}^C = q_{z=0}^C + Q^L ,\end{equation}

which is consistent with our prediction (4.18).

Figure 10. Heat balance variation with respect to $R$ (a), $\lambda _\mu$ (b) and $K$ (c). Each term has been normalised by the heat flux in at the bottom boundary $q^c_{z=0}=\theta$. In (a) $K = 1000, \lambda _\mu = 10^6$; in (b) $K = 1000, R = 10^6$; in (c) $R = 5\times 10^5, \lambda _\mu = 2 \times 10^5$.

To interpret the results we first recall from (4.8) and (4.5) that the latent heat $Q^L \sim R ^{1/2} [ \phi ^s u_z^s ]_{z=1}$. As discussed above, the solid velocity increases with $R$ as $u_z^s \sim R^{1/2}$ at low $R$ and tends to an $R$-independent limit as $R\to \infty$ while the average solid volume fraction varies as $\phi ^s \sim R^{-1/2}$ as $R\to \infty$. Thus, the dependence of $Q^L$ on $R$ cancels out for sufficiently large $R$. As a result, $Q^L$ increases with $R$ but becomes independent of $R$ as $R\to \infty$ (figure 10a). As more latent heat is absorbed for melting, less heat flows out at the top of the layer. The temperature gradient becomes less steep (less negative), and the temperature difference between top and bottom decreases (see figure 13 in Appendix E). Magnitudes of heating due to heat pipe $Q^P$ and friction $Q^F$ increase with $R$; both terms are proportional to $R^{1/2}$ (recall (4.8)), while for low $R$ the solid flux increases $j_z \sim R^{1/2}$. At large $R, j_z \sim R^{-1/2}$ and thus the heating terms $Q^P$ and $Q^F$ become independent of $R$. The viscous heating term $Q^V$ increases initially with $R$ due to increased gradients of the phase velocity near the bottom boundary (recall figure 2c,d). However, since $\phi ^s \sim R^{-1/2}$ as $R\to \infty$, the viscous heating due to the solid phase decreases for large $R$, and thus the overall viscous heating $Q^V$ decreases with $R$.

The effect of increasing $\lambda _\mu$ has an opposite effect to that of increasing $R$. As $\lambda _\mu$ is increased the solid-phase velocity is decreased, tending to $u_z^s \sim \lambda ^{-1}$ as $\lambda _\mu \to \infty$. Thus, the latent heat $Q^L$ decreases as $\lambda _\mu$ increases (figure 10b). Smaller $Q^L$ implies less heat lost to melting, which results in a steeper (more negative) temperature gradient and larger temperature difference across the layer. Similarly, $Q^P$ and $Q^F$ decrease with $\lambda _\mu$ owing to the decrease in the solid flux ($\, j_z \sim \lambda _\mu ^{-1}$). Viscous heating $Q^V$ also decreases for very large $\lambda _\mu$, as increased viscous forces smooth out all sharp velocity gradients, thus reducing the amount of viscous heating.

As $K$ is increased, magnitudes of velocities of both phases as well as the solid flux are decreased ($u_z^s \sim u_z^l \sim j_z \sim K^{-1}$ as $K\to \infty$). Thus, all contributions to the heat balance decrease as $K$ is increased (figure 10c).

5.6. Effects of compaction viscosity

Thus far we considered solutions with constant (shear) viscosity. As described in § 2.6, a compaction viscosity $\mathcal {C}$ also arises in two-phase systems due to the resistance of solid to compaction. Here, we briefly examine the effect of $\mathcal {C}$ on the solutions and the limiting behaviour at large $\lambda _\mu$. Little is known about the functional form of $\mathcal {C}$ for the iron in the F-layer. At high solid fraction a common functional form is $\mathcal {C} \sim 1/\phi ^l$ (e.g. Bercovici et al. 2001). However, for two-phase materials that deform via diffusion creep, the form $\mathcal {C} \sim - \log \phi ^l$ may be more appropriate (e.g. Rudge 2018; Katz et al. 2022). Thus, we opt for a simple logarithmic form:

(5.3)\begin{equation} \mathcal{C}(\phi^l) =- C_0\log \phi^l, \end{equation}

where $C_0$ is a dimensionless parameter measuring the relative magnitude of compaction viscosity to shear viscosity. Compaction viscosity $\mathcal {C}(\phi ^l)$ is zero in the absence of solid phase ($\phi ^l = 1$), and has a logarithmic singularity at vanishing porosity ($\phi ^l \to 0$).

Compaction viscosity $\mathcal {C}(\phi ^l)$ enters the system in three places: in the solid momentum equation, in the form of an additional stress term; in the liquidus relation, where it affects the melting temperature; and in the temperature equation in the form of an additional heating term. As explained in § 4.3, we expect the latter two effects to be negligible. In the solid momentum equations (4.1c), compaction viscosity combines with the shear viscosity to make the effective viscosity of the solid $\lambda _\mu ( 4/3 - C_0 \log \phi ^l)$. Thus, increasing $C_0$ increases the solid viscosity, which reduces $u_z^s$ (making it less negative) and increases $\phi ^s$ and $u_z^l$. Indeed, this can be observed in figure 11, which shows profiles of the freezing rate $\varGamma ^s$, the solid volume fraction $\phi ^s$ and vertical velocities of solid and liquid phases $u_z^s, u_z^l$, for increasing values of $C_0$. As $C_0$ is increased, solid phase becomes more compact (porosity $\phi ^l$ decreases), and liquid phase is expelled at faster velocities. Solid-phase velocities decrease, leading to decreased solid flux and concurrent reduction in melting required to balance that mass flux.

Figure 11. Profiles of freezing rate $\varGamma ^s$ (a), solid volume fraction $\phi ^s$ (b), solid-phase velocity $u_z^s$ (c) and liquid-phase velocity $u_z^l$ (d) for increasing values of the compaction parameter $C_0$; other parameters are fixed at $R=10^6, \lambda _\mu = 10^4, K = 100$.

Figure 12 shows the variation of average solid velocity, average liquid velocity and average solid volume fraction versus $\lambda _\mu$ for different values of $C_0$. For greater $C_0$, effective viscosity is greater and so the large-$\lambda _\mu$ regime is entered at lower $\lambda _\mu$. The scaling behaviour $u_z^s\sim \phi ^l\sim \lambda _\mu ^{-1}$ as $\lambda _\mu \to \infty$ seems to persist in the presence of porosity-dependent viscosity.

Figure 12. Variation of the average solid-phase velocity $K|\langle u_z^s \rangle | / (\lambda _\rho -1)$ (a), the average liquid-phase velocity $K\langle u_z^l \rangle / (\lambda _\rho -1)$ (b) and the average liquid volume fraction $\langle \phi ^l \rangle$ (c) with respect to the viscosity ratio $\lambda _\mu$ for different values of compaction parameter $C_0$; other parameters are fixed at $R=10^6, K = 100$. In (a,b) the phase velocities have been scaled by the estimated terminal velocity $(\lambda _\rho - 1)/K$.

6. Discussion and conclusions

In this paper we have developed and analysed a model of two-phase slurry dynamics, motivated by recent studies of the F-layer at the base of Earth's outer core. The governing equations we solve are an extension of the two-phase equations derived by Roberts & Loper (1987) to include pressure disequilibrium between the two phases (Bercovici et al. 2001; Ricard 2007). The equations contain a more general formulation of the fluid dynamics compared with those employed in previous studies of solid–liquid layers in planetary cores (e.g. Rückriemen et al. 2015; Davies & Pommier 2018; Wong et al. 2018; Bercovici & Mulyukova 2021) and can be applied to a wide range of two-phase systems including ‘iron snow’ and ‘flotation’ layers in the cores of small terrestrial bodies.

We have studied 1-D time-independent solutions of the two-phase equations in order to establish systematic behaviour in the parameter regime $R \gg 1, \lambda _\mu \gg 1$ relevant to planetary cores. A key benefit of the 1-D system is that the numerous control parameters can be varied over many orders of magnitude. Considering, separately, the variation in the nature of the solutions as $R\to \infty$, or $\lambda _\mu \to \infty$, reveals two distinct regimes of limiting behaviour. In the limit of $R\to \infty$ the solutions are characterised by small values of solid fraction and liquid-phase velocity $\phi ^s\sim u_z^l \sim R^{-1/2}$, while the solid-phase velocity tends to the Darcy velocity $u_z^s \sim -(\lambda _\rho -1)/K$. On the other hand, in the limit of $\lambda _\mu \to \infty$, the solutions are characterised by small values of liquid fraction and solid-phase velocity $\phi ^l\sim u_z^s \sim \lambda _\mu ^{-1}$, while the liquid-phase velocity tends to the Darcy velocity $u_z^l \sim (\lambda _\rho -1)/K$. The transition between the two regimes of limiting behaviour in $R$–$\lambda _\mu$ space occurs roughly when $R$ and $\lambda _\mu$ are comparable in magnitude. The developed scalings are in excellent agreement with our numerical solutions of the full slurry equations.

In the parameter regime of planetary cores where $R > \lambda _\mu \gg 1$ the average solid fraction in the layer is small. Near the bottom boundary all solutions exhibit a thin region with $\phi ^s \to 1$; however, in the parameter regime of interest ($R > \lambda _\mu \gg 1$), the relative thickness of this region to the depth of the entire layer is of the order of 1 % (see for example the case in figure 2b with $R = 10^8$). Extrapolation of our results to F-layer conditions (figure 9b) indicates that a pure iron slurry F-layer would contain a mean solid fraction of at most 5 %. This result is consistent with Gubbins et al. (2008), who used a purely thermodynamical line of reasoning. Such low solid fractions are incompatible with seismically inferred densities, supporting the view that the F-layer in Earth's core is not a pure iron slurry.

There are several physical effects that may arise in Earth's F-layer that are not captured within our current model. Earth's liquid core contains around 10 wt% lighter elements such as O, Si and S and previous models suggest that compositional variations are likely important for explaining the stable stratification of the F-layer (Gubbins et al. 2008; Wong et al. 2021). Generalising our pure iron model to multiple chemical constituents is relatively straightforward (Keller & Suckale 2019), but comes at the expense of adding more input parameters to the 14 already present and also introducing new physical effects. In a multicomponent slurry, additional buoyancy forces will arise from differences in specific densities and concentrations of individual constituents, which can facilitate compositional convection. Furthermore, when density depends on two properties (here, heat and light element concentration) that diffuse at different rates, double diffusive processes can occur. Nevertheless, we have shown that the pure slurry displays a wide range of dynamical behaviour, including solutions with low solid fraction across much of the layer that will serve as an important reference point for future studies of the multicomponent system.

Our model does not include the effects of rotation or magnetic field, both of which are known to be important for the dynamo process in planetary cores. However, in the stratified F-layer, their dynamical significance is much less clear. Incorporating rotation into the governing equations is conceptually trivial, it merely involving inclusion of the Coriolis effect in the momentum equations of each phase. Incorporating magnetic fields into the two-phase equations in not a trivial task, but it is possible (Bercovici & Mulyukova 2020). However, in both cases the analysis will inevitably become more complicated as it calls for fully three-dimensional treatment, which is beyond the scope of the present study.

When developing the governing equations we introduced a number of simplifications with the aim of reducing mathematical complexity while still retaining the key physical processes. The fast melting approximation considerably simplifies the thermodynamics; it assumes that both phases share the same temperature and constrains the system to remain in phase equilibrium, where the temperature follows the liquidus. Thermal equilibrium is expected to be very rapid, $O(1)$ s, which is much faster than the time scales associated with the flow velocities we have obtained and so equality of solid and liquid temperatures appears to be a reasonable assumption. The effect of departures from phase equilibrium is much harder to assess. Walker et al. (2021) modelled the non-equilibrium effects of nucleation and crystal grow of falling particles in the F-layer and their preliminary results indicate that the layer would not reach phase equilibrium. However, the physics of crystal nucleation and growth in planetary cores is complex and is only now starting to be investigated in detail (Huguet et al. 2018; Davies, Pozzo & Alfè 2019; Wilson et al. 2021, 2023) so it is presently unclear whether the departures from phase equilibrium would have a significant effect on the fluid dynamics. Departures from phase equilibrium could be incorporated into the governing equations using macroscopic parametrisations of the microscopic nucleation and growth processes (Loper 1992), but at the expense of significantly complicating the theory. At present we believe the fast-melting limit is a reasonable compromise.

A significant challenge in modelling two-phase liquid–solid flows lies in capturing the rheological behaviour of the solid phase. Real liquid–solid mixtures are known to exhibit complex rheological behaviour which is a function of the solid volume fraction and can be either Newtonian or non-Newtonian in nature (Guazzelli & Pouliquen 2018). In this paper we adopted a simple Newtonian rheology with uniform shear viscosity for both phases. This was partly to understand the fundamental properties of the system and partly because there is little information available about the rheology of solid iron at core conditions. Clearly more work is desirable to constrain the functional form of phase viscosities in a slurry. Nevertheless, in the slurry limit (small $\phi ^s$), modifications to the viscosity are of the order $\phi ^s$ (e.g. Guazzelli & Pouliquen 2018) and thus are not expected to cause drastic changes in the nature of solutions.

In deriving the two-phase equations it was necessary to postulate functional forms relating velocities and pressures of the two phases to the interfacial velocity and pressure. The associated coefficients $\gamma$ and $\omega$ were prescribed based on simple physical considerations, whereas in reality they should be determined from experimental and/or computational studies of the microscale behaviour of the two-phase system. Clearly more work is desirable in order to better constrain these quantities, which can then be applied in a straightforward manner to the theory we have derived. However, since the associated terms are neglected in the reduced model we have developed, which represents the gross behaviour of our numerical solutions very well, we believe that the specific functional forms of these coefficients do not critically impact our results.

The 1-D basic state is the first step in mathematical analysis of fluid dynamical systems, it being the underlying equilibrium (stable or not) on top of which three-dimensional perturbations evolve. In single-phase fluid flow problems, the underlying basic state tends to be trivial, with fluid at rest in hydrostatic balance. On the other hand, in a two-phase pure slurry, even the 1-D steady state is interesting because the two phases move relative to one another. This 1-D state is part of any time-evolving three-dimensional solution, and it is thus integral to understand its nature. We believe that the detailed analysis performed in this study constitutes a successful first step towards these longer-term goals and paves the way for analytical considerations of the more complete model of the slurry.