2. Laboratory experiments and key observations

Previously described laboratory experiments provide a motivation and context for testing the theory developed here. We focus on experiments conducted on partially molten rock, typically synthesized from mixtures of $\sim$95 % olivine grains and $\sim$5 % mid-ocean ridge basalt, sometimes with a small percentage of chromite (e.g. Holtzman et al. Reference Holtzman, Groebner, Zimmerman, Ginsberg and Kohlstedt2003; King et al. Reference King, Zimmerman and Kohlstedt2010; Qi et al. Reference Qi, Kohlstedt, Katz and Takei2015). The olivine grains are polydisperse, typically with a mean diameter of $\sim$10 $\mathrm {\mu }$m. Samples are hydrostatically hot pressed to remove gas-filled bubbles prior to deformation. After hot pressing, they have a nominally uniform melt fraction, $\phi _0$.

The experiments are conducted in a gas-medium, triaxial-deformation apparatus (Paterson Reference Paterson1990) with a confining pressure of 300 MPa and temperatures of ${\sim }1225\,^\circ$C. The samples are jacketed to separate them from the confining gas. Under these conditions, the basalt is molten and the olivine (and chromite) grains are solid. Torsional deformation is imposed on the sample, although some experiments have also been conducted in direct shear (Holtzman et al. Reference Holtzman, Groebner, Zimmerman, Ginsberg and Kohlstedt2003; Holtzman & Kohlstedt Reference Holtzman and Kohlstedt2007). The distribution of porosity within the sample is not measured in situ. Rather, the experiment is quenched, sectioned and imaged at high resolution to calculate the porosity field.

The essential characteristics of these torsional experiments are outlined in figure 1(a). Cylindrical samples with height $H$ and outer radius $R$ are deformed by a circular plate that turns with angular velocity $\dot {\varOmega }\hat {\boldsymbol {z}}$ about the axis of the cylinder. At low strains, when the sample remains nominally uniform, the imposed twist induces an azimuthal velocity field $U(r,z)\hat {\boldsymbol {\varphi }}$ that is axisymmetric. The velocity component $U$ increases linearly in the $\hat {\boldsymbol {z}}$ and $\hat {\boldsymbol {r}}$ directions. On cylindrical surfaces, the deformation is approximately that of simple shear; the magnitude of the shear strain (and its rate) increase from zero at the twist axis to a maximum value at the outer boundary of the sample.

Figure 1. Experimental configuration and representative results. (a) Schematic diagram of a deforming experimental sample and the emergent patterns of melt segregation. Experiments are conducted at high confining pressure and high temperature. After achieving a specified twist, the sample is quenched, sectioned and polished to reveal the distribution of melt (solidified to glass) and crystalline, granular solid. (b) A tangential section showing high-porosity bands (black) at low angle to the shear plane ($\phi _0=0.04, \gamma =1.5$; King et al. Reference King, Zimmerman and Kohlstedt2010). (c) A transverse section showing radially inward migration of melt ($\phi _0=0.10, \gamma =5.0$; Qi et al. Reference Qi, Kohlstedt, Katz and Takei2015). Cracks visible in (b) and (c) are a consequence of the rapid quench and decompression after deformation.

The outer boundary of the sample is sealed in an impermeable, nickel jacket. The radial normal stress at the jacket is maintained constant by the confining gas pressure. At the temperatures of the experiments, the viscosity of nickel is greater than that of the basaltic liquid and less than the granular olivine aggregate. This viscosity contrast enables the jacket to shear with negligible resistance, but discourages its intrusion into the pore space of the sample. Hence, its effect on the sample falls somewhere between two limiting cases. In one limit, the jacket inhibits all radial flow at the boundary by isolating the volume of the sample. In the other limit, it transmits the full confining pressure into the pore space between olivine grains. There is empirical evidence that the reality is closer to the first of these limits, but the details have not been measured or quantified.

Two critical observations have arisen from torsion experiments on partially molten rocks. The first is the emergence of high-porosity sheets separated by compacted, low-porosity lenses after a shear strain of $\gamma \sim 1$ (Holtzman et al. Reference Holtzman, Groebner, Zimmerman, Ginsberg and Kohlstedt2003). These sheets are usually measured in cross-section, as in figure 1(b), and hence, referred to as bands. They have a characteristic spacing and form at 15–20$^\circ$ to the shear plane (Holtzman & Kohlstedt Reference Holtzman and Kohlstedt2007). This angle is similar to that of Riedel shear zones (Dresen Reference Dresen1991), but significantly lower than what would be expected if the bands were normal to the direction of maximum tension (45$^\circ$). Furthermore, individual bands are embedded in a nominally simple-shear flow and, hence, with time, they are rotated to higher angles. However, despite this necessary rotation, the band-angle distribution remains roughly unchanged with increasing strain (King et al. Reference King, Zimmerman and Kohlstedt2010).

The second critical observation is that with progressive twist, liquid melt segregates from the solid grains and migrates toward the centre of the cylinder (Qi & Kohlstedt Reference Qi and Kohlstedt2018). This migration leads to azimuthally averaged porosity, measured over the transverse section shown in figure 1(c), that decreases with radius. Experiments to increasing values of total twist (reported as shear strain $\gamma$ at the outer radius) exhibit greater segregation and a steeper radial porosity gradient (Qi et al. Reference Qi, Kohlstedt, Katz and Takei2015).

The present study aims to explain these observations in terms of the physics of dense granular suspensions.

3. Rheological model

Our rheological model of a two-phase aggregate, comprising a contiguous matrix of solid grains and its melt-saturated, permeable pore space, is based on the poro-viscous theory derived by McKenzie (Reference McKenzie1984) and reviewed by Katz (Reference Katz2022). The melt is present with volume fraction $\phi$ (the porosity) that varies in space and time. This variation is accommodated by (de)compaction of the solid matrix, but both phases are incompressible. To incorporate dilatancy effects, we take inspiration from theories of suspensions (Brady & Morris Reference Brady and Morris1997; Fang et al. Reference Fang, Mammoli, Brady, Ingber, Mondy and Graham2002; Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018) and append a term that hypothetically quantifies the normal stresses generated by grain–grain interactions during shearing flow. The constitutive law for the effective stress is then

(3.1)\begin{equation} \boldsymbol{\sigma}^{eff} = \zeta_\phi\mathcal{C}\boldsymbol{I} + 2\eta_\phi\dot{\boldsymbol{\varepsilon}} - D_\phi\boldsymbol{\varLambda}\dot{\varepsilon}_{II}, \end{equation}

where $\zeta _\phi$, $\eta _\phi$ and $D_\phi$ are dynamic viscosities for isotropic, deviatoric and dilatational deformation, respectively, $\boldsymbol {I}$ is the identity tensor, and where

(3.2a–c)\begin{equation} \mathcal{C} \equiv \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}^s,\quad \dot{\boldsymbol{\varepsilon}} \equiv \frac{1}{2}\left[\boldsymbol{\nabla}\boldsymbol{v}^s + (\boldsymbol{\nabla}\boldsymbol{v}^s)^T - \frac{2}{3}\mathcal{C}\boldsymbol{I}\right]\!,\quad \boldsymbol{\varLambda}\equiv \left(\begin{array}{@{}ccc@{}} 1 & 0 & 0 \\ 0 & \varLambda_\perp & 0 \\ 0 & 0 & \varLambda_\times \end{array}\right)\!, \end{equation}

are the decompaction rate, deviatoric strain-rate tensor and particle-stress anisotropy tensor, respectively. We have introduced $\boldsymbol {v}^s$, the solid velocity field, and $\dot {\varepsilon }_{II} \equiv \sqrt {\dot {\boldsymbol {\varepsilon }}:\dot {\boldsymbol {\varepsilon }}/2}$ is the second invariant of the deviatoric strain-rate tensor. Deboeuf et al. (Reference Deboeuf, Gauthier, Martin, Yurkovetsky and Morris2009) provides theoretical context, insightful commentary and empirical justification for the dilatancy term in (3.1).

The particle-stress anisotropy tensor $\boldsymbol {\varLambda }$ is used to model the normal stresses generated by a particle-laden flow that is locally approximated as simple shear (Guazzelli & Morris Reference Guazzelli and Morris2011; Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018). It is written with reference to a coordinate system aligned with the simple shear. The $\varLambda _{11}$ direction is taken to be the direction of flow (indicated by $\parallel$); the $\varLambda _{22}$ direction is normal to the shear plane (hence, we denote it $\varLambda _\perp$); the $\varLambda _{33}$ direction is the direction of the vorticity vector (and, hence, denoted $\varLambda _\times$). The $\varLambda _{11}$ entry is factored out and lumped with $D_\phi$. Therefore, $\varLambda _\perp$ and $\varLambda _\times$ are dimensionless particle-normal-stress ratios. Previous work has shown that the values of these parameters may be constrained by comparison of model predictions with carefully designed experiments (Morris & Boulay Reference Morris and Boulay1999; Fang et al. Reference Fang, Mammoli, Brady, Ingber, Mondy and Graham2002; Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018). In the flow geometries considered below (Cartesian or cylindrical), the particle-stress anisotropy tensor can be straightforwardly aligned with the experimental deformation geometry; in general, it must be aligned with respect to the principal axes of the flow (Miller, Singh & Morris Reference Miller, Singh and Morris2009).

The isotropic part of the effective stress is where dilatancy modifies the physics. We see this by taking $-1/3$ the trace of the effective stress tensor in (3.1),

(3.3)\begin{equation} D_\phi\dot{\varepsilon}_{II}\, \text{tr}\left(\boldsymbol{\varLambda}\right)/3 = (1-\phi)(P^s-P^\ell) + \zeta_\phi\mathcal{C}, \end{equation}

where $P^j = -\text {tr}\left (\boldsymbol {\sigma }^j\right )/3$ is the pressure of phase $j$. This equation states that the shear-strain rate has two possible consequences for isotropic deformation. If $\zeta _\phi =0$, there is no viscous resistance to compaction and shear generates a positive effective pressure. This is equivalent to suspension theory (e.g. Deboeuf et al. Reference Deboeuf, Gauthier, Martin, Yurkovetsky and Morris2009). If, in contrast, there is zero effective pressure (${\rm \Delta} P=0$), then shear causes dilatation. This has a parallel in soil mechanics, where the dilatancy angle $\psi$ gives the kinematic relationship between shear strain and dilatation (e.g. Oda & Iwashita Reference Oda and Iwashita2020). In the present context, we can compute the dilatancy angle as $\tan \psi \equiv D_\phi \text {tr}\left (\boldsymbol {\varLambda }\right )/3\zeta _\phi$. The more general case, of interest here, is where both the effective pressure and the compaction viscosity are non-zero.

To complete the rheological model, we require expressions for the dependency of the three viscosities on melt fraction $\phi$. Empirical constraints and theoretical models of the shear viscosity $\eta _\phi$ and compaction viscosity $\zeta _\phi$ are summarised by Katz et al. (Reference Katz, Rees Jones, Rudge and Keller2022). Shear viscosity has been measured over a range of melt fractions; Kelemen et al. (Reference Kelemen, Hirth, Shimizu, Spiegelman and Dick1997) showed that it is well described by an exponential decrease with liquid fraction $\phi$. Theory for Coble creep, where compaction is accommodated by diffusion of grain mass along grain boundaries and through the melt-filled pores, indicates that the compaction viscosity is a multiple of $\sim$5/3 larger than the shear viscosity (Takei & Holtzman Reference Takei and Holtzman2009; Rudge Reference Rudge2018). There are no empirical measurements of the dilitation viscosity $D_\phi$ of partially molten rock, nor are there are microstructural models. Experiments on particle suspensions by Deboeuf et al. (Reference Deboeuf, Gauthier, Martin, Yurkovetsky and Morris2009) show an exponential weakening of particle normal stress with liquid fraction. On this basis, and for simplicity in the absence of further information, we take $D_\phi$ to be an unknown multiple of $\eta _\phi$. Hence, the viscosities are given by

(3.4a–c)\begin{equation} \eta_\phi = \eta_0\,\text{e}^{-\lambda(\phi-\phi_0)},\quad \zeta_\phi = 5\eta_\phi/3,\quad D_\phi=D_0\eta_\phi, \end{equation}

where $\eta _0$ is a reference value of shear viscosity at reference melt fraction $\phi _0$, $\lambda \approx 27$ is the porosity-weakening factor and $D_0$ is an unknown, dimensionless constant.

We obtain a constraint on $D_0$ by requiring positive entropy production under any combination of shear and isotropic deformation. The dissipation-rate density arising from (3.1) is

(3.5)\begin{equation} \varPsi = \zeta_\phi\mathcal{C}^2 + 4\eta_\phi\dot{\varepsilon}_{II}^2 - D_\phi\dot{\varepsilon}_{II}(\dot{\varepsilon}_\parallel+ \varLambda_\perp\dot{\varepsilon}_\perp+ \varLambda_\times\dot{\varepsilon}_\times), \end{equation}

where $\dot {\varepsilon }_\parallel, \dot {\varepsilon }_\perp, \dot {\varepsilon }_\times$ are the normal components of the strain-rate tensor in a coordinate system aligned with simple shear. Assuming isotropic dilatancy $\boldsymbol {\varLambda } = \boldsymbol {I}$, (3.5) becomes $\varPsi = \zeta _\phi \mathcal {C}^2 + 4\eta _\phi \dot {\varepsilon }_{II}^2 - D_\phi \dot {\varepsilon }_{II}\mathcal {C}$, and into this we substitute the viscosities of (3.4a–c). We find that $\varPsi$ is positive definite if $0 \le D_0 < 4\sqrt {5/3} \approx 5$ and, therefore, limit consideration to values of $D_0$ within this range.

Finally, in combining our rheological model with conservation equations governing the flow, we consider the granular physics discussed by Kamrin & Koval (Reference Kamrin and Koval2012), which adopts a model for emulsions by Goyon et al. (Reference Goyon, Colin, Ovarlez, Ajdari and Bocquet2008). They show that macroscopic, irreversible shear is accommodated by grain-rearrangement events at the microscopic scale. In partially molten rock, geometric compatibility of the grain packing dictates that grains cannot rotate freely; their rotations must be compatible with those of neighbouring grains (Rudge Reference Rudge2021). Hence, deformation is necessarily dispersed by grain–grain interaction during rearrangement events. This non-local interaction means that the viscosity at a point in the medium is influenced by the viscosity at points within a distance $\xi$, known as the cooperativity length. Kamrin & Koval (Reference Kamrin and Koval2012) express this interaction in terms of a non-local fluidity – the inverse of the non-local shear viscosity $\tilde {\eta }_\phi$. We rewrite their fluidity equation in terms of a non-local viscosity,

(3.6)\begin{equation} \tilde{\eta}_\phi^{-1} = \eta_\phi^{-1} + \xi^2\nabla^2 (\tilde{\eta}_\phi^{-1}). \end{equation}

Evidently, if $\xi =0$ then the non-local viscosity reduces to $\eta _\phi$. For $\xi >0$, this equation imposes a minimum scale of viscosity variation. Goyon et al. (Reference Goyon, Colin, Ovarlez, Ajdari and Bocquet2008) measure $\xi$ as a function of $1-\phi$ and find that it increases to about 5 times the grain diameter at a solid fraction of 85 %. We shall see below in § 4.2 that $\xi >0$ serves to regularise the spectrum of instability growth.

4. Analysis

To explore the consequences of the hypothesised rheological model, we adopt the formulation of mass and momentum conservation for a partially molten rock deforming at zero Reynolds number (McKenzie Reference McKenzie1984). We make a Boussinesq approximation, taking density as constant for both phases, assume zero mass transfer between phases and neglect gravitational body forces on the basis that they are much weaker than the shear tractions imposed in experiments. With these assumptions, the phase densities vanish from the equations. The coupled system of conservation equations becomes

(4.1a)$$\begin{gather} \mathcal{C} = \boldsymbol{\nabla}\boldsymbol{\cdot} M_\phi \boldsymbol{\nabla} P^\ell, \end{gather}$$

(4.1b)$$\begin{gather}\boldsymbol{\nabla} P^\ell = \boldsymbol{\nabla}\boldsymbol{\cdot}2\tilde{\eta}_\phi\dot{\boldsymbol{\varepsilon}} + \boldsymbol{\nabla}\tilde{\zeta}_\phi\mathcal{C} - \boldsymbol{\nabla}\boldsymbol{\cdot}\tilde{D}_\phi\boldsymbol{\varLambda}\dot{\varepsilon}_{II}, \end{gather}$$

(4.1c)$$\begin{gather}\frac{\text{D}_{s}{\phi}}{\text{D}{t}} = (1-\phi)\mathcal{C}. \end{gather}$$

The first, known as the compaction equation, is obtained from Darcy's law by eliminating the liquid velocity using the two-phase continuity equation. It includes the fluid mobility $M_\phi = M_0 (\phi /\phi _0)^n$, which represents the ratio of the porosity-dependent permeability and the constant liquid viscosity. The second equation is a statement of force balance in the two-phase aggregate. The third equation is mass conservation for the solid phase (porosity is transported by the solid velocity). These equations are standard (Katz Reference Katz2022), except for two modifications. The first modification is use of the non-local viscosity $\tilde {\eta }_\phi$ in (4.1b), which couples it to (3.6) governing the non-local viscosity. For consistency, we use the non-local viscosity in (3.4a–c) to compute non-local compaction $\tilde {\zeta }_\phi$ and dilitation $\tilde {D}_\phi$ viscosities. The second modification is the last term on the right-hand side of (4.1b), which captures the hypothesised dilatancy effects. The classical model is recovered for ${\xi,D_0\to 0}$.

In the subsections below, we investigate the consequences of these two modifications. We do so in the context of torsional deformation and boundary conditions that mimic the laboratory experiments described in § 2.

4.1. Radial segregation in parallel-plate torsion

Torsional flow embeds simple shear into a cylindrical geometry with the potential for hoop stress. The experiments described in § 2 demonstrate that parallel-plate torsional flow drives solid radially outward and liquid radially inward. This phenomenon is consistent with the behaviour of dense suspensions undergoing parallel-plate torsional flow (Merhi et al. Reference Merhi, Lemaire, Bossis and Moukalled2005) but in contrast to the Poiseuille flow of Appendix D. We consider cone-and-plate torsional flow (where the plates are not parallel) in Appendix C.

To understand the radially outward transport of solid grains in terms of dilatancy and particle-stress anisotropy, we work in a cylindrical geometry with coordinates ($r,\varphi,z$), as shown in figure 1(a). We consider a cylinder of partially molten rock with outer radius $R$, azimuthal symmetry in $\varphi$ and, instantaneously at $t=0$, with uniform porosity $\phi _0$. At this instant, the solid flow is assumed to have zero $\hat {\boldsymbol {z}}$ component, a fixed azimuthal component and an unknown radial component. This flow is described by

(4.2a–c)\begin{equation} \boldsymbol{v}^s = V\hat{\boldsymbol{r}} + \frac{\dot{\varOmega}\,r z}{H}\hat{\boldsymbol{\varphi}}, \quad\mathcal{C}=\frac{1}{r}\frac{\partial{}}{\partial{r}}rV,\quad\dot{\boldsymbol{\varepsilon}} = \left(\begin{array}{@{}ccc@{}} \displaystyle \dfrac{\partial{V}}{\partial{r}} - \dfrac{\mathcal{C}}{3} & 0 & 0 \\ 0 & \dfrac{V}{r} - \dfrac{\mathcal{C}}{3} & \dfrac{\dot{\varOmega r}}{2H} \\ 0 & \dfrac{\dot{\varOmega r}}{2H} & -\dfrac{\mathcal{C}}{3} \end{array}\right)\!, \end{equation}

where $V(r)$ is the unknown radial component of the solid velocity field, $\dot {\varOmega }$ is the constant twist rate and $H$ is the uniform gap between the parallel plates. We linearise the strain-rate intensity under the assumption that dilatancy is driven by the forced shear such that

(4.3)\begin{equation} \dot{\varepsilon}_{II}\sim\dot{\varOmega}r/2H. \end{equation}

This choice eliminates a feedback whereby the anisotropic part of the dilatancy drives additional dilatancy. While it may be physically reasonable, it will reduce the predicted dilatancy at a given value of $D_0$ relative to the case where the feedback is included.

We use $\dot {\varepsilon }_{II}$ in the radial component of force-balance equation (4.1b) to write

(4.4)\begin{equation} \frac{\partial{P^\ell}}{\partial{r}} = 3\eta_0\frac{\partial{}}{\partial{r}}\frac{1}{r}\frac{\partial{}}{\partial{r}}rV - \frac{D_0\dot{\varOmega}R^2}{6H}(2\varLambda_\times-1). \end{equation}

We then combine this with the compaction equation (4.1a) to eliminate the liquid pressure and integrate once. Rescaling $r$ with the outer radius $R$ and $V$ with the characteristic scale

(4.5)\begin{equation} [V] = \frac{D_0\dot{\varOmega}R^2}{6H}, \end{equation}

we obtain the dimensionless equation

(4.6)\begin{equation} \frac{\partial{}}{\partial{r}}\frac{1}{r}\frac{\partial{}}{\partial{r}}rV - \frac{V}{\mathcal{R}^2} = 2\varLambda_\times-1. \end{equation}

Here we have introduced

(4.7)\begin{equation} \mathcal{R} \equiv \frac{\sqrt{3\eta_0M_0}}{R}, \end{equation}

the ratio of the compaction length to the outer radius. The compaction length is an emergent length scale over which perturbations to the solid–liquid pressure difference are relaxed by decompaction (McKenzie Reference McKenzie1984; Spiegelman Reference Spiegelman1993; Katz Reference Katz2022).

The normal-stress difference on the right-hand side of (4.6) arises from the particle-stress anisotropy tensor $\boldsymbol {\varLambda }$, with the coordinates aligned such that the flow direction is $\hat {\boldsymbol {\varphi }}$ and the vorticity direction is $\hat {\boldsymbol {r}}$. The boundary condition at the centre of the cylinder is $V(0)=0$. With this constraint, (4.6) admits the solution

(4.8)\begin{equation} V(r) = {\rm \pi}\mathcal{R}^2(\tfrac{1}{2}-\varLambda_\times)[A I_1(r/\mathcal{R}) - L_1(r/\mathcal{R})], \end{equation}

where $I_n(z)$ is the modified Bessel function of the first kind, $L_n(z)$ is the modified Struve function and $A$ is a constant to be determined by matching the boundary condition at the dimensionless outer radius $r=1$. Two end-member cases can be considered for this outer boundary condition.

4.1.1. Outer boundary condition: no normal flow

In this case, a rigid outer cylinder requires that at $r=R$, the radial component of velocity is $V(R)=0$. Then the analytical solution to dimensionless equation (4.6) is

(4.9)\begin{equation} V(r) = {\rm \pi}\mathcal{R}^2\left(\frac{1}{2} -\varLambda_\times\right)\left[ \frac{L_1(1/\mathcal{R})}{I_1(1/\mathcal{R})}I_1(r/\mathcal{R})-L_1(r/\mathcal{R})\right]\!. \end{equation}

This result demonstrates that the sign of $V(r)$ is determined by the size of $\varLambda _\times$. In figure 2(a) we have chosen $\varLambda _\times = 0.45$ such that $V(r)>0$. This choice is qualitatively consistent with experimental results (§ 2) where the solid is observed to move outward with progressive twist. Evidently, for the outer boundary condition $V(1)=0$, any choice of $\varLambda _\times$ satisfying

(4.10)\begin{equation} \varLambda_\times < 1/2 \end{equation}

is also qualitatively consistent. For this range of $\varLambda _\times$, the hoop stress generated by dilatancy in the flow direction is stronger than the dilatant normal stress in the radial (vorticity) direction. As noted by Takei & Katz (Reference Takei and Katz2013), a compressive hoop stress drives solid radially outward.

Figure 2. Parallel-plate torsion flow at $t=0$. Non-dimensional solutions of (4.6) with uniform $\eta _\phi =\eta _0$. Panels (a,b) have the outer boundary condition $V(R)=0$; panels (c,d) have the zero-effective-stress outer boundary condition given in (4.12). (a) Analytical solutions (4.9) with the outer boundary condition $V(R)=0$ and with $\varLambda _\times =0.45$. In the limit of $\mathcal {R} \gg 1$, the dimensionless solution is asymptotic to $V(r)\sim (2\varLambda _\times - 1) (r^2-r)/3$. In the other limit, $\mathcal {R}\ll 1$, the matched asymptotic solution is $V(r)\sim \mathcal {R}^2(2\varLambda _\times - 1)[-1+\exp (-r/\mathcal {R}) + \exp (-(1-r)/\mathcal {R})]$. (b) Decompaction rate with $\varLambda _\times =0.45$. (c) Analytical solution (4.13) with outer boundary condition (4.12) and with $\mathcal {R}=0.3$. (d) Decompaction rate with $\mathcal {R}=0.3$.

4.1.2. Outer boundary condition: no normal effective stress

Alternatively, we can consider the case where the partially molten cylinder is surrounded by an inviscid fluid, held at a dimensional confining pressure $P_c$. This pressure must be balanced by the phase-averaged traction at the boundary and, therefore, $\hat {\boldsymbol {r}} \boldsymbol {\cdot } \bar {\boldsymbol {\sigma }}\boldsymbol {\cdot } \hat {\boldsymbol {r}} = -P_c$. Assuming that the liquid pressure is continuous at $r=R$, we obtain the boundary condition

(4.11)\begin{equation} \hat{\boldsymbol{r}}\boldsymbol{\cdot} \boldsymbol{\sigma}^{eff} \boldsymbol{\cdot} \hat{\boldsymbol{r}} = 0 \quad\text{ at } r=R. \end{equation}

Expanding this condition using (3.1) and the approximate second invariant (4.3), then non-dimensionalising $r$ with $R$ and $V$ with $[V]$, we obtain the dimensionless boundary condition

(4.12)\begin{equation} \frac{V}{3} + \frac{\partial{V}}{\partial{r}} = \varLambda_\times \quad\text{ at } r=1. \end{equation}

This condition yields a different value of $A$ and the general solution (4.8) becomes

(4.13)\begin{align} V(r) &= {\rm \pi}\mathcal{R}^2\left(\frac{1}{2}-\varLambda_\times\right)\notag\\ &\quad\times \left[\frac{3L_0\left(\dfrac{1}{\mathcal{R}}\right) - 2\mathcal{R} L_1 \left(\dfrac{1}{\mathcal{R}}\right)+\dfrac{3}{{\rm \pi}\mathcal{R}} \dfrac{\varLambda_\times}{1/2-\varLambda_\times}}{\displaystyle 3I_0 \left(\dfrac{1}{\mathcal{R}}\right) - 2\mathcal{R} I_1\left(\dfrac{1}{\mathcal{R}}\right)}I_1 \left(\frac{r}{\mathcal{R}}\right)-L_1\left(\frac{r}{\mathcal{R}}\right)\right]. \end{align}

This function is plotted in figure 2(c,d); the curves are computed with $\mathcal {R}=0.3$ and values of $\varLambda _\times$ that span $1/2$. The radial component of solid velocity $V$ is generally positive, indicating outward solid flow and decompaction across all radii. This outward flow is again driven by the compressive hoop stress. Distinct from the rigid outer boundary condition, however, condition (4.12) allows the solid to move outward at the outer boundary. Hence, in this case, torsion with any $\varLambda _\times$ causes the solid cylinder to expand radially, imbibing liquid across the outer boundary.

The pattern of flow in figure 2(c) is slightly different for $\varLambda _\times =0.75$, where there is a region of $V<0$ at inner radii. The inner part of this region is associated with compaction $\mathcal {C}<0$, as shown by the solid curve in (d). The driving force is again dilatancy, but in this case with $\varLambda _\times >1/2$, the radial dilatant normal stress plays a significant role. As is the case for Poiseuille flow in Appendix D, faster shear at larger radii drives solid inward. But with zero radial effective stress at the outer boundary, dilatancy also drives outward solid flow, radial expansion of the cylinder and radial imbibition of liquid.

4.1.3. Outward force on pistons due to dilatancy

The results depicted in figure 2 for both outer boundary conditions are valid instantaneously at $t=0$, when all properties are uniform with radius. At this initial instant, we compute an axial force outward on the plates, parallel to $\hat {\boldsymbol {z}}$, that arises from dilatancy. This calculation provides a prediction to be compared with laboratory measurements. Details of the calculation are in Appendix A. The main result is that the axial force is dominated by the direct effect of dilatancy in the flow-perpendicular direction, and hence, scales as

(4.14)\begin{equation} {\rm \Delta} F \approx \mathcal{T} D_0 \varLambda_\perp{/} 3R, \end{equation}

where $\mathcal {T}$ is the torque that causes a twist rate of $\dot \varOmega$ at $t=0$. Here ${\rm \Delta} F$ is the outward force in excess of that due to the confining pressure surrounding the sample. Using typical laboratory values for $\mathcal {T}$ and $R$ in (4.14) (Appendix A), we find that the excess force is of the order of 10 % of the force due to the typical confining pressure.

In detail, the excess force ${\rm \Delta} F$ can deviate from the simple prediction of (4.14). Figure 8 in Appendix A plots this deviation for both outer-boundary-condition cases over a range of $\mathcal {R}$ and $\varLambda _\times$. When the outer boundary is closed to solid flow ($V(R)=0$), the excess force is close to the simple scaling above; when the outer boundary has zero effective stress (4.12), dilatancy in the radial direction leads to a net decompaction of the sample and a reduction in the excess axial force by approximately one half.

4.1.4. Finite time and steady state

The analysis of parallel-plate torsion to this point has considered the instantaneous problem at $t=0$, when the domain is uniform in porosity. The instantaneous flow requires that this uniform state is subsequently lost by radial segregation of solid and liquid (and, as shown empirically and below in § 4.2, by a banding instability). Appendix B derives the system of equations, simplified from (4.1), that governs the finite-time evolution of the radial distribution of porosity. The non-uniform porosity $\phi (r)$ leads to a radial dependence of mobility $M_\phi$ and aggregate viscosity $\tilde {\eta }_\phi$.

A series of time-dependent numerical solutions are plotted in figure 3(a), coloured according to the shear strain $\gamma$ at the outer radius of the domain $R$. Details of the numerical method are in Appendix B; the code is available in an online repository (Katz, Rudge & Hansen Reference Katz, Rudge and Hansen2023). This calculation uses $\mathcal {R}=0.3$ and $\varLambda _\times =0.4$ with boundary condition $V(R)=0$. At the smallest finite strains, the porosity distribution has the shape of the $t=0$ solution for $\mathcal {C}(r)$, shown in figure 2(b). With increasing strain (time), the porosity contrast between the centre and outer radius increases. However, the evolution slows and ceases as the porosity distribution approaches a steady state.

Figure 3. Parallel-plate torsion at $t\ge 0$ and $t\to \infty$. Coloured curves show the time-dependent, numerical solution to the system (B3) for porosity $\phi (r,t)/\phi _0$. Black curves show the analytical, steady-state solution (4.16) for $\xi =0$. Both panels use empirically motivated values $\lambda =27$ and $\phi _0=0.07$. (a) Solutions with $\varLambda _\times =0.4$, $\mathcal {R}=0.3$ and $\xi /R=0.03$ at various outer-radius strains $\gamma (R)$. (b) Steady solutions for four values of $\varLambda _\times$.

In the steady state and with $\xi =0$, the radial component of the solid velocity is zero and radial force-balance equation (4.1b) reduces to $\hat {\boldsymbol {r}}\boldsymbol {\cdot } (\boldsymbol {\nabla }\boldsymbol {\cdot } D_\phi \boldsymbol {\varLambda }\dot {\varepsilon }_{II}) = 0$ or, after expanding and rearranging,

(4.15)\begin{equation} D_\phi = \varLambda_\times\left(2D_\phi + rD_\phi'\frac{\partial{\phi}}{\partial{r}}\right)\!, \end{equation}

where $D_\phi '$ is the derivative of the dilatation viscosity with respect to its argument, $\phi$. In solutions to this equation, a steady state is reached when the force of the hoop stress (left-hand side) balances the force of the radial normal stresses (right-hand side). The compressive hoop force, associated with a coefficient of unity in $\boldsymbol {\varLambda }$, is due to azimuthal dilatancy that pushes solid radially outward. The radial normal force, associated with the coefficient $\varLambda _\times$ in $\boldsymbol {\varLambda }$, has two causes: first, the gradient in radial dilatancy due to the torsional shear ($\dot {\varepsilon }_{II}\propto r$), and second, the gradient in radial dilatancy due to the steady-state radial gradient in porosity.

Laboratory experiments that impose torsional deformation, discussed in § 2, have an azimuthally averaged porosity that decreases with radius. On the basis of the predicted compaction rate at $t=0$, shown in figure 9, we can infer that this porosity structure is consistent with the no-normal-flow boundary condition and $\varLambda _\times < 1/2$. It is unclear whether these experiments approach a steady-state radial porosity profile, or would do so at larger strains. If a steady state can be achieved, then (4.15) requires that $D_\phi '<0$ independent of the specific form of $D_\phi$; in other words, it requires that the dilatancy stress at fixed strain rate decreases as porosity increases.

In (3.4a–c) we specified that $D_\phi$ has the form $D_\phi =D_0\eta _\phi$. Given the exponential dependence of $\eta _\phi$ on $\phi$, it follows that $D_\phi ' = -\lambda D _\phi$. With this we can solve (4.15) to give

(4.16)\begin{equation} \phi(r,t\to\infty) = \phi_0 - \frac{1/2-\varLambda_\times}{\varLambda_\times\lambda/2}\left(\ln r + \frac{1}{2}\right)\!, \end{equation}

where $r$ has been non-dimensionalised with the outer radius $R$ and we have used global conservation of liquid mass to determine the constant of integration (see Appendix B for details). This function is plotted as black curves in figure 3 for various values of $\varLambda _\times$; in § 5 we compare it with measurements from laboratory experiments. The logarithmic singularity in (4.16) for $r \to 0$ is removed by non-local viscosity when cooperativity length $\xi >0$. This is evident in figure 3(a) by comparison between the numerical solution at late time (red curve; $\xi =0.1\delta$) and the steady solution (black curve; $\xi =0$).

4.2. Simple shear between parallel plates

Here, following the analysis of Spiegelman (Reference Spiegelman2003), we investigate the stability of a two-dimensional simple-shear flow with initial melt fraction $\phi _0 + \phi _1(\boldsymbol {x},t)$, where $\vert \phi _1\vert \ll \phi _0$ is a perturbation. A schematic diagram is shown in figure 4(a). The coordinate system is oriented such that $\hat {\boldsymbol {x}}$ is in the flow direction and $\hat {\boldsymbol {y}}$ is in the direction perpendicular to the shear plane. We assume invariance in the $\hat {\boldsymbol {z}}$ (vorticity) direction and take the $x$–$y$ plane to be infinite; hence, there is no need to impose boundary conditions. The procedure is a standard linearised stability analysis, detailed in Katz (Reference Katz2022, Chapter 7) and sketched in the next paragraph. Alisic et al. (Reference Alisic, Rhebergen, Rudge, Katz and Wells2016) provides a three-dimensional analysis for torsion in cylindrical coordinates, but this adds mathematical complexity without additional physical insight.

Figure 4. Growth rate of sinusoidal perturbations under a simple-shear flow from (4.17). (a) Schematic diagram showing a finite region of the infinite domain. The grey scale shows the perturbed porosity field.(b) Growth rate as a function of wavenumber $k$ for $D_0=0$ at $\theta =45^\circ$. Circles represent the growth rate computed at $k^*=\epsilon ^{-1/2}$. (c) Growth-rate angular factor as a function of $\theta$ with $\varLambda _\perp = 1$, and values of $D_0$ given in the legend. (d) Growth-rate angular factor as a function of $\theta$ with $D_0=2$, and values of $\varLambda _\perp$ given in the legend.

We use (4.1b) to eliminate the pressure gradient from (4.1a) and obtain an equation governing the irrotational part of the velocity field. Then we take the curl of (4.1b) to obtain an equation governing the solenoidal part. These are coupled to (3.6) and (4.1c) for viscosity and solid mass. We expand variables into a steady, background state and a time-dependent perturbation that is arbitrarily small at $t=0$. The perturbations are assumed to be proportional to $\exp [{\rm i}\boldsymbol {k}(t)\boldsymbol {\cdot } \boldsymbol {x} + {s}(t)]$, where $\boldsymbol {k}(t)$ is a time-dependent wave vector that changes direction and magnitude with the background flow. After linearising the governing equations, we solve the leading-order balance for the base state. This is a simple-shear flow with velocity gradient $\dot \gamma$ and zero compaction rate. Using the base-state solution, we solve the perturbation equations to obtain the dimensionless growth rate $\dot {s}$ as a function of dimensionless wave-vector magnitude $k$ and wavefront angle to the shear plane $\theta \equiv \tan ^{-1}k_x/k_y$. In the present case, we obtain

(4.17)\begin{equation} \dot{s} = (1-\phi_0)\frac{\lambda}{3}\frac{k^2}{(1+k^2)(1+\epsilon^2 k^2)}\left[\sin 2\theta - \frac{D_0}{2}\frac{(\sin^2\theta + \varLambda_\perp\cos^2\theta)\sin^2 2\theta}{1-\dfrac{D_0}{4}(1-\varLambda_\perp)\sin 4\theta}\right]\!. \end{equation}

In this equation, $\dot {s}$ has been made dimensionless by scaling with the background rate of shear $\dot {\gamma }$. Wavenumber $k$ has been made dimensionless by scaling with the compaction length $\delta \equiv \sqrt {3 M_0 \eta _0}$. We have introduced the ratio $\epsilon \equiv \xi /\delta$ representing the dimensionless cooperativity scale. Localisation phenomena in partially molten rock can emerge at scales smaller than the compaction length. In this case, localisation occurs because positive perturbations have lower viscosity, and hence, decompact under resolved tension (Stevenson Reference Stevenson1989; Spiegelman Reference Spiegelman2003). We refer to these perturbations as ‘bands’, making explicit reference to the high-porosity bands seen in experimental cross-sections. Below we discuss the dependence of band growth rate on wavenumber, angle and physical parameters.

Figure 4(b) shows the wavenumber dependence of the growth rate for several values of $\epsilon$ (assuming optimal orientation, $\theta = \theta ^*$). The curve for $\epsilon =0$ is the case with zero cooperativity of the viscosity field. It shows the classical result that all wavelengths smaller than the compaction length ($k\gg 1$) grow equally fast (Stevenson Reference Stevenson1989). The use of the non-local viscosity with finite $\epsilon$ regularises the spectrum, imposing a short-wavelength cutoff at a dimensional wavelength $\sim \xi$, the cooperativity scale of the non-local viscosity. Growth-rate curves in figure 4(b) have a maximum at a dimensional wavenumber

(4.18)\begin{equation} k^* = 1/\sqrt{\xi\delta}. \end{equation}

Figure 4(c) shows the growth rate as a function of the angle $\theta$ between wavefronts and the shear plane. Four curves show different values of the dilatancy viscosity prefactor $D_0 \ge 0$ (assuming $\boldsymbol {\varLambda } = \boldsymbol {I}$, i.e. isotropic dilatancy). The curve for $D_0=0$ corresponds to the case with no dilatancy, as studied by Spiegelman (Reference Spiegelman2003). This case has positive growth rates between zero and 90$^\circ$, a range over which bands are subject to tension, and negative rates for angles greater than 90$^\circ$, which are subject to compression. The maximum growth rate occurs at $\theta ^*=45^\circ$, where band wavefronts are perpendicular to the principal tension axis.

For larger $D_0$ in figure 4(c), dilatancy leads to peak growth rate at low and high angles. In particular, the two maxima of growth rate occur at angles $\theta ^*_{1,2}$ that vary with $D_0>1$,

(4.19a,b)\begin{equation} {\theta}^*_1 = \arcsin(1/D_0)/2,\quad {\theta}^*_2 = {\rm \pi}/2 - \arcsin(1/D_0)/2. \end{equation}

A growth-rate peak at $\theta ^*=15^\circ$, which is roughly that observed in experiments, corresponds to $D_0=2$.

Dilatancy has two competing effects that combine to produce this spectral shift with increasing $D_0$. First, due to $D_\phi '<0$, the background simple-shear flow causes a dilatancy perturbation that is exactly anti-phase with the porosity perturbation. This causes perturbations to decay at a rate that is independent of band angle $\theta$. Second, the porosity perturbations create variations in shear viscosity, which in turn create perturbations in the rate of shear strain. These drive variations in dilatancy that are exactly in-phase with variations in porosity; they hence contribute to perturbation growth. Critically, however, the growth rate associated with this second mechanism depends on band angle $\theta$. Shear localises when bands are at low or high angle to the shear plane and, therefore, this effect is proportional to $\cos ^22\theta$. The combination of the two contributions of isotropic dilatancy is negative, overall, and proportional to $\sin ^22\theta$.

Dilatancy in partially molten rock may be anisotropic, however. Experiments on dense granular suspensions, reviewed by Guazzelli & Pouliquen (Reference Guazzelli and Pouliquen2018), have obtained inconclusive estimates of $\varLambda _\perp$, with some reporting $\varLambda _\perp$ increasing from unity with solid fraction, some reporting the opposite, and others reporting $\varLambda _\perp$ within error of unity for all solid fractions. In figure 4(d) we compare growth-rate curves as a function of band angle for $\varLambda _\perp = 0.8, 1, 1.2$ for $D_0=2$. The differences in the growth-rate peaks are subtle – likely indistinguishable on the basis of measured band-angle histograms from experiments.

5. Comparison with laboratory data

Published results from laboratory experiments (§ 2) provide an opportunity to test our theory. We begin by considering the best-established outcome of experiments, the localisation of melt into high-porosity bands. In particular, we first consider the angle that the bands make to the shear plane. In figure 5 blue points with error bars indicate the mean and standard deviation of band angles from experiments quenched at shear strains between about 1 and 4. Each panel presents the same data set. The points form a coherent array with mean angles in the range 15–20$^\circ$, except at strains greater than about 3, at which the points spread out into a range of 10–25$^\circ$. The dashed lines, also identical in each panel, indicate trajectories that band angles would follow if rotated passively in the simple-shear flow (Katz et al. Reference Katz, Spiegelman and Holtzman2006).

Figure 5. Angle spectra of porosity-band amplitude as a function of shear strain $\dot {\gamma }t$. The data points are the same in each panel. They record the mean angle from band-angle histograms of individual, published experiments (see legend); error bars are one standard deviation of the histogram. Solid lines are contours of the band amplitude $\exp [s(t)]$, normalised over angles at each increment of strain. Dashed lines are passive advection trajectories (see text). Amplitude is computed by quadrature of the growth rate ${\dot {s}(t)}$ from (4.17) with $\varLambda _\perp = 1$ and dimensionless $k(t=0) = k^* = \epsilon ^{-1/2}$. Each panel has a different magnitude of dilatancy: (a) $D_0=0$;(b) $D_0=2$; (c) $D_0=3$.

Comparison of the array of blue points with the adjacent dashed lines clearly demonstrates that the evolution cannot be characterised as passive rotation of an initial set of bands. The maintenance of low angles over large strains has been attributed to successive generations of low-angle bands that draw melt from (and hence replace) previous generations as they undergo rotation to angles unfavourable for growth (Holtzman, Kohlstedt & Morgan Reference Holtzman, Kohlstedt and Morgan2005; Katz et al. Reference Katz, Spiegelman and Holtzman2006). This process occurs at a finite perturbation amplitude and, strictly speaking, should not be described by solutions of linearised governing equations. However, if the angular spectrum of growth rate (i.e. figure 4c) remains approximately independent of strain, then a forward integration of $\dot {s}$ with respect to time (strain) along the trajectories of passive rotation might approximate the evolving angular spectrum of amplitude. In this context, normalisation of the amplitude spectrum at each increment of strain might qualitatively represent melt redistribution.

The contours of this finite-strain, normalised perturbation amplitude spectrum are plotted in figure 5 for three different values of $D_0$ (each with $\boldsymbol {\varLambda } = \boldsymbol {I}$). In (a), for $D_0=0$, we see that without dilatancy, peak predicted amplitudes are far from observations. In (b), for $D_0=2$, and in (c), for $D_0=3$, we see that peak predicted amplitudes occur at angles close to measured values. The $D_0=2$ case is in better alignment at lower strains where nonlinear effects might be less important, and hence, we take this value of the dilatancy prefactor to be most appropriate in the context of the present model assumptions.

The linearised theory for porosity-band growth also provides a prediction of the wavelength with the largest growth rate. Again, this prediction is strictly valid only near the onset of instability when the porosity contrast remains small. Laboratory experiments that produce bands provide the opportunity to measure mean band width and spacing; summing these provides an estimate of the band wavelength. However, as with band angles, these measurements are made in a nonlinear regime when the porosity contrast is large. So a comparison with theory is not strictly valid. However, as with band angles, if the growth-rate spectrum remains roughly independent of strain, then a correspondence between theory and experiment might be expected even at larger strains. In that case, we would expect the observed wavelength to be proportional to $1/k^* = \sqrt {\xi \delta }$, where $\xi$ is the cooperativity length scale and $\delta$ is the compaction length.

Figure 6 suggests that this correspondence may hold. The blue symbols represent the wavelength of bands from experiments as a function of the geometric mean of the empirically known grain diameter $d$ and compaction length $\delta$. The cooperativity scale $\xi$ is understood to be proportional to grain diameter (Henann & Kamrin Reference Henann and Kamrin2013) and, hence, $1/k^* \propto \sqrt {d\delta }$. Although there is considerable uncertainty on the experimental estimates (particularly $\delta$), a linear trend is compatible with the data.

Figure 6. Wavelength of porosity bands in laboratory experiments (see legend) plotted against the geometric mean of the grain size $d$ and the compaction length. The data-source publications provide mean estimates of band spacing, band width, grain size and compaction length. Band wavelength is calculated as the sum of mean band spacing and mean band width. Errors on measurements are propagated to give the error bars. The dashed line is a fit to the data respecting uncertainties on both axes (York et al. Reference York, Evensen, Martinez and De Basabe Delgado2004; Wiens Reference Wiens2023).

Finally, we evaluate model predictions of the radial distribution of porosity in torsion experiments quenched at different strains. The experiments are a subset of those published by Qi et al. (Reference Qi, Kohlstedt, Katz and Takei2015) and Qi & Kohlstedt (Reference Qi and Kohlstedt2018); we select only those in which the liquid phase is basaltic melt. We re-analyse high-resolution binary images of transverse sections, following the authors’ published protocol, but averaging azimuthally over fewer, wider rings. Recalculated porosities are presented as a function of normalised sample radius in figure 7. The data are coloured by the magnitude of shear strain at the outer radius, which ranges from 0 to $\sim$14. Three experiments with strains between 5 and 6 are averaged to produce one radial series. Error bars represent the standard deviation of the averaged and normalised porosity at $t=0$, at which the porosity should be uniform.

Figure 7. Radial distribution of porosity $\phi (r)$ normalised by the initial porosity $\phi _0$ in experiments and theory. Symbols represent porosity from laboratory experiments obtained by reprocessing high-resolution scans of transverse sections (Qi et al. Reference Qi, Kohlstedt, Katz and Takei2015; Qi & Kohlstedt Reference Qi and Kohlstedt2018); values are averages over rings of equal radial span. Error bars show the standard deviation of porosity in the undeformed ($\gamma =0$) experiment. Colours represent the shear strain at the outer radius $\gamma (R)$. The black, dotted line represents the steady-state solution (4.15) with $\varLambda _\times =0.4$, $\phi _0 = 0.04$, $\lambda =27$. Dashed curves are numerical solutions to the system (B3) with the same parameters as for the steady curve and also $D_0=2$, $\mathcal {R}=0.3$. The numerical solution is plotted at finite values of $\gamma (R)$, as given by their colour. Appendix B gives details of the numerical method.

A qualitative conclusion can be immediately drawn by examining the data. The boundary condition imposing zero effective stress at $r=R$ is not consistent with the experiments. This is because of the pattern of decompaction shown in figure 2(d), which has the most rapid decompaction at the outer boundary. In contrast, the empirical data uniformly exhibit reduced porosity due to compaction there. Hence, we proceed using only the $V(R)=0$ condition of no radial flow at the outer boundary.

Model predictions are overlayed onto the laboratory data in figure 7. The black dotted line is the steady-state solution from (4.16). This is computed with $\varLambda _\times = 0.4$, a value chosen to give an approximate match with the highest-strain experiment. (Note, however, that we have no evidence that the empirical porosity distribution at $\gamma \approx 14$ is in steady state.) The dashed curves are numerical solutions of the time-dependent model (Appendix B), plotted at values of outer-radius strain indicated by the colour of the curve. The shape of the model curves is in qualitative agreement with the data trends, within error. However, model porosity evolves more rapidly as a function of strain than the porosity in experiments.

There are two main difficulties in interpreting this mismatch in terms of parameter values or deficiencies of the theory. The first is that we do not know whether the porosity distributions from experiments shown in figure 7 are approaching a steady state, as predicted by the theory. The uncertainties in the porosity and the coarse sampling in total strain make any inferences speculative. Second, supposing that the experiments are evolving toward a steady state, we do not have a reliable estimate of the porosity distribution in that state. If we had constraints on that distribution (and knowledge of $\phi _0$ and $\lambda$), we could use (4.16) to infer $\varLambda _\times$.

Advancing speculatively, we assume that the experiments do approach a steady state. We note that smaller $\varLambda _\times$ corresponds to larger $\text {d}\phi /\text {d} r$ at steady state (figure 3b). Assuming that our estimates for $\phi _0\approx 0.04$ and $\lambda \approx 27$ are sufficiently accurate, the data in figure 7 are indicative of $\varLambda _\times \lesssim 0.4$. We can estimate the time scale of porosity adjustment to steady state in the numerical solutions shown in figure 7. Referring to the porosity evolution (4.1c), we approximate $\text {D}_s\phi /\text {D}t$ by ${\rm \Delta} \phi /\tau$, where $\tau$ is the time scale over which porosity changes by ${\rm \Delta} \phi$ from its initial value $\phi _0$ to its steady value at some given radius. According to (4.1c), this change is driven by decompaction at a rate $(1-\phi )\mathcal {C}$; we approximate this rate as $\breve {\mathcal {C}}$, which we take to be the decompaction rate at $r=0, t=0$ (c.f. figure 2b). This rate is obtained by calculating $\mathcal {C}(r)$ at $t=0$ from the analytical solution (4.9) for radial velocity, re-dimensionalising and evaluating at $r=0$ with $\mathcal {R} \sim 1$. We then form the outer-radius strain at time $\tau$ as $\gamma (\tau ) = \tau /\dot {\gamma } = {\rm \Delta} \phi /\dot {\gamma }\breve {\mathcal {C}}$, where $\dot {\gamma } = \dot {\varOmega }R/H$ is the outer-radius strain rate. This obtains

(5.1)\begin{equation} \gamma(\tau) \sim \frac{18}{D_0\varLambda_\times\lambda}, \end{equation}

where we have used (4.16) to evaluate ${\rm \Delta} \phi$ at $r = 1/\text {e}^2$. For the parameters used in figure 7, this gives an outer-radius strain of $\gamma \sim 1$ at which the simulated porosity has evolved to within a factor of $1/\text {e}$ of its steady value. This estimate is comparable to the numerical solution of figure 7, but it is smaller than the empirical time scale by a factor of 10. We cannot bring these time scales into agreement by changing $D_0$ because this is constrained by the angle of porosity bands (figure 5), and we have already assumed that we know $\lambda$. Reducing $\varLambda _\times$ thus appears to be an option. According to figure 3(b), reducing $\varLambda _\times$ by a factor of 2 predicts a steady-state, radial porosity gradient much larger than observed at the largest empirically attained strain $\gamma$. Experiments to larger $\gamma$ are needed to test this prediction. Alternatively, the discrepancy in time scales may be a consequence of nonlinear interaction of radial segregation with emergence of high-porosity layers, which is not captured in our models.