Granular dilatancy and non-local fluidity of partially molten rock

Abstract

Partially molten rock is a densely packed, melt-saturated, granular medium, but it has seldom been considered in these terms. In this paper we extend the continuum theory of partially molten rock to incorporate the physics of granular media. Our formulation includes dilatancy in a viscous constitutive law and introduces a non-local fluidity. We analyse the resulting poro-viscous–granular theory in terms of two modes of liquid–solid segregation that are observed in published torsion experiments: localisation of liquid into high-porosity sheets and radially inward liquid flow. We show that the newly incorporated granular physics brings the theory into agreement with experiments. We discuss these results in the context of grain-scale physics across the nominal jamming fraction at the high homologous temperatures relevant in geological systems.

1. Introduction

Partially molten rock is a physical system that is central to many geological and planetary processes. It is a densely packed, melt-saturated, granular medium, but it has seldom been considered in these terms. Continuum models of partially molten rock treat the solid and liquid phases as interpenetrating fluids in a poro-viscous, zero-Reynolds-number theory (e.g. McKenzie 1984; Fowler 1990). In such models, effects arising from the discrete grains are neglected, except insofar as they affect the creep viscosity. For example, during Coble creep, the melt phase provides a fast pathway for mass diffusion around grains (e.g. Takei & Holtzman 2009; Rudge 2018). Deformation experiments on partially molten rock are typically parameterised in terms of an isotropic flow law with a weakening factor that depends on the volume fraction of melt within pores (e.g. Kohlstedt & Zimmerman 1996; Kelemen et al. 1997). These physics were reviewed by Katz et al. (2022).

It is less commonly noted, however, that deformation of partially molten rock inevitably includes a component of sliding along grain boundaries (e.g. Hansen, Zimmerman & Kohlstedt 2011; Rudge 2021). The granular origins and importance of such sliding in partially molten rock were recognised by Paterson (1995) and elaborated by Paterson (2001). In these works, the geometric grain-compatibility problem arising from grain-boundary sliding is assumed to be entirely resolved by shape change of the grains, which occurs by diffusion or lattice dislocations (Langdon 2006). A third possibility is noted by Paterson (2001) but then neglected: that incompatibility is resolved by relative motion of undeforming grains in a granular flow. This mechanism is fundamental in the physics of athermal granular media (e.g. Forterre & Pouliquen 2008); it gives rise to dilatancy and non-local granular fluidity.

In this paper we extend the continuum theory of partially molten rock to incorporate the physics of a granular medium. As a hypothesis for the essential granular physics, we adapt and include theory for dilatancy and non-local fluidity. We test this hypothesis by modelling published laboratory experiments in which partially molten rock is subjected to torsional deformation. The deformation drives liquid–solid segregation and yields robust patterns of melt localisation. We show that the inclusion of granular physics brings model predictions into agreement with laboratory data.

The laboratory experiments, detailed in King, Zimmerman & Kohlstedt (2010) and reviewed in § 2 below, are conducted on synthetic rocks comprising solid olivine grains and liquid basaltic melt. Hot-pressed, nominally uniform, cylindrical samples of this aggregate are sheared in a torsion apparatus at high temperature and confining pressure. During shear, two modes of liquid–solid segregation occur simultaneously. The first is a pattern-forming localisation of the liquid into high-porosity sheets that form at 15–20$^\circ$ to the shear plane (Holtzman et al. 2003). The sheets are typically measured in their cross-section, where they appear as bands with a characteristic spacing. The second is a radially inward porous flow of liquid, accommodated by a radially outward flow of solid (Qi et al. 2015).

A satisfactory physical understanding of these flow phenomena has been elusive, though much has been learned through theoretical analysis. Localisation of the liquid phase into sheets at 45$^\circ$ to the shear plane was predicted by Stevenson (1989) and Spiegelman (2003) to be a consequence of a porosity-weakening viscosity of the solid aggregate. Katz, Spiegelman & Holtzman (2006) and Rudge & Bercovici (2015) showed that if this viscosity is (effectively) non-Newtonian with a power-law exponent of $\sim$6, the angle is reduced to match the observations. However, this exponent was measured by King et al. (2010) to be ${\sim }1.5\pm 0.3$ at 95 % confidence – almost Newtonian. Furthermore, these isotropic theories cannot explain the radial melt segregation in torsion experiments. In contrast, a theory of anisotropic Coble creep with Newtonian viscosity can explain the radial segregation (Qi et al. 2015). This theory was derived by Takei & Holtzman (2009) from grain-scale considerations of anisotropic solid contiguity under deviatoric stress. It predicts that viscous resistance to deformation is reduced in the direction of minimum contiguity. It also predicts the emergence of porosity bands and, if the contiguity tensor aligns with the principal stress directions, that the bands grow fastest at low angles, consistent with experiments (Takei & Katz 2013). However, in laboratory experiments that produce bands, the grain-scale contiguity is misaligned by about 15$^\circ$ (Qi et al. 2015; Qi & Kohlstedt 2018), which corresponds to a theoretical prediction of high-angle porosity bands (this discrepancy is resolved by better measurement of contiguity, according to Seltzer et al. 2023). Nonetheless, viscous anisotropy theory gives rise to an effective dilatancy that we discuss in § 6.1, below.

Another challenge is to explain the characteristic wavelength of the high-porosity bands observed in experiments. All of the theories noted above lack mode selection; instead, they predict the rate of band growth to plateau at decreasing wavelength. Several studies have invoked processes driven by interfacial energy to regularise the growth-rate spectrum. Bercovici & Rudge (2016) incorporated capillary effects in a diffuse-interface approximation of a sharp porosity interface (Sun & Beckermann 2004) – however, sharp interfaces emerge in experiments only long after the onset of instability. Takei & Hier-Majumder (2009) and King, Hier-Majumder & Kohlstedt (2011) hypothesised that variation of surface tension drives dissolution/precipitation reactions. When coupled with chemical diffusion in the melt phase, these reactions damp instability growth at small wavelengths.

The theory of dense granular suspensions, as reviewed by Guazzelli & Pouliquen (2018), holds promise in providing a simple and unified explanation for all of these observed patterns. A central feature is the anisotropic compressive stress between solid particles caused by shearing flow (Bagnold 1954). The coupling between shear and compression is the consequence of the microphysical interaction of suspended particles (Brady & Morris 1997). This behaviour is demonstrated empirically in various studies, but Deboeuf et al. (2009) provide a particularly fascinating example and discussion. If the suspended solid phase is not rigidly confined, it can undergo a net dilatation due to shear (Reynolds 1885; Boyer, Guazzelli & Pouliquen 2011).

In a suspension contained within a constant volume, net dilatancy is prohibited but the solid fraction can vary internally. Besseling et al. (2010) shows that shearing, dense suspensions are susceptible to a banding instability; this instability is modelled in terms of a suspension viscosity that increases with solid fraction. Their results run parallel to the theory for band emergence in partially molten rock (Stevenson 1989) except that in a suspension, the growth rate of bands also depends on the dilatancy. Moreover, Morris & Boulay (1999) shows that for suspension flows in cylindrical geometry (i.e. pipe flow, parallel-plate or cone-and-plate torsion), radial segregation of liquid and solid phases is predicted, consistent with experiments. Again, there is a parallel with results for partially molten rock in torsion (Qi et al. 2015; Qi & Kohlstedt 2018) and pipe-flow (Quintanilla-Terminel et al. 2019) configurations. In all of these flowing suspensions, the dilatancy stress plays a central role.

Another aspect of granular physics that may be relevant here is non-local fluidity (inverse viscosity). This concept was developed in the context of emulsions (e.g. Goyon et al. 2008; Bocquet, Colin & Ajdari 2009) and adapted to granular suspensions (Kamrin & Koval 2012). The theory states that the flow response to stress at a point in the granular medium is sensitive to the fluidity in a neighbourhood around that point. This neighbourhood has a typical size, $\xi$, of order $10$ times the grain size and decreasing with the square root of shear stress. Henann & Kamrin (2013) demonstrate that simulated shear zones, forced by a spatial discontinuity in boundary velocity, are regularised by non-local fluidity to a width that is consistent with experiments. Hence, non-local fluidity appears promising in regularising the growth spectrum of shear bands in experiments on partially molten rock.

The fundamentally granular nature of partially molten rock and the relevant predictions from theories of dense granular suspensions motivate the present work. Our aims are to develop a theory for partially molten rock that incorporates granular dilatancy and non-local fluidity, and to compare predictions of that theory to the results of laboratory experiments. We note that in doing so, we are applying granular physics at solid fractions above what is typically considered the jamming fraction, at which the solid phase becomes immobile. However, in crystalline materials at high homologous temperatures, grain boundaries behave as a viscous fluid that allows grains to slide past each other (Ashby 1972). In this context, the solid phase can still be mobilized if sufficient shear stress is applied (Heussinger & Barrat 2009). Moreover, in-situ observations of polycrystalline aggregates deforming at high solid fraction show a clear link between grain-boundary sliding and dilatancy (Walte, Bons & Passchier 2005; Kareh et al. 2017). Hence, we assert that although the grains of partially molten rocks are not rigid, they nonetheless undergo grain-boundary sliding that is associated with a compressive intergranular stress and may lead to dilatancy.

Mechanical decreases in solid fraction, including by dilatancy, are here referred to as decompaction. The poro-viscous theory of partially molten rock relates the decompaction rate to the pressure difference between the liquid and solid phases in a viscous constitutive law (McKenzie 1984). This approach differs from suspension theory, in which the solid phase exerts zero resistance to changes in solid fraction (Guazzelli & Morris 2011). It also differs from theories for dry granular media, where the solid fraction is a decreasing function of the shear-strain rate (Forterre & Pouliquen 2008), and from soil mechanics, where the solid fraction is predicted to evolve toward a critical state as a function of the total strain (Oda & Iwashita 2020). However, it seems that these isotropic dynamics may be incompletely understood. For example, Kabla & Senden (2009) found empirical evidence that dilatancy and shear-independent compaction compete in the evolution of solid fraction. By combining poro-viscous decompaction with dilatancy stress, our theory may provide new insight in this regard.

Previous authors have incorporated granular dilatancy into discussions and models of geological materials, going back at least to Mead (1925). It has been invoked in crystal-rich deforming magma (e.g. Smith 1997; Petford, Koenders & Clemens 2020), in lower-crustal shear zones (Menegon et al. 2015) and in gouge-filled fault zones (e.g. Marone, Raleigh & Scholz 1990; Segall et al. 2010). Dilatation has been considered in competition with compaction (Paterson 2001; Niemeijer & Spiers 2007) and as a microphysical mechanism responsible for rate-and-state friction (Chen & Spiers 2016). It may play a role in regulating glacial sliding (Warburton, Hewitt & Neufeld 2023) and in a range of geomorphological processes (Jerolmack & Daniels 2019). Dilatation is associated with Riedel shear zones (e.g. Dresen 1991; Bedford & Faulkner 2021), which appear at the same angle as bands in partially molten rock. It might be expected that partially molten rock shares certain behaviours with other granular, geological materials. In the present work we find that incorporation of granular dilatancy and non-local fluidity brings predictions of a poro-viscous compaction theory into quantitative agreement with experimental results.

The paper is organised as follows. In § 2 we review torsion experiments on partially molten rocks and highlight their key results. We present our rheological model in § 3. Then, in § 4 we provide the governing equations and analyse them in terms of radial segregation and band formation. This analysis is followed by quantitative comparison with experiments in § 5 and a discussion in § 6.

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. 2003; King et al. 2010; Qi et al. 2015). 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 1990) 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. 2003; Holtzman & Kohlstedt 2007). 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. 2010). (c) A transverse section showing radially inward migration of melt ($\phi _0=0.10, \gamma =5.0$; Qi et al. 2015). 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. 2003). 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 2007). This angle is similar to that of Riedel shear zones (Dresen 1991), 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. 2010).

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 2018). 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. 2015).

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 (1984) and reviewed by Katz (2022). 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 1997; Fang et al. 2002; Guazzelli & Pouliquen 2018) 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. (2009) 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 2011; Guazzelli & Pouliquen 2018). 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 1999; Fang et al. 2002; Guazzelli & Pouliquen 2018). 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 2009).

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. 2009). 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 2020). 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. (2022). Shear viscosity has been measured over a range of melt fractions; Kelemen et al. (1997) 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 2009; Rudge 2018). 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. (2009) 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 (2012), which adopts a model for emulsions by Goyon et al. (2008). 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 2021). 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 (2012) 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. (2008) 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 1984). 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 2022), 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. 2005) 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 1984; Spiegelman 1993; Katz 2022).

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 (2013), 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 2023). 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 (2003), 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 (2022, Chapter 7) and sketched in the next paragraph. Alisic et al. (2016) 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 1989; Spiegelman 2003). 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 1989). 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 (2003). 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 (2018), 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. 2006).

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 2005; Katz et al. 2006). 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 2013) 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. 2004; Wiens 2023).

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. (2015) and Qi & Kohlstedt (2018); 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. 2015; Qi & Kohlstedt 2018); 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.

6. Discussion

We hypothesised that granular physics (i.e. dilatancy and non-local fluidity) shapes the patterns that emerge when partially molten rock is deformed in laboratory experiments. Our model predictions, based on a rheological formulation combining theories for poro-viscous compaction and dense granular suspensions, can be made quantitatively consistent with most aspects of the empirical data. This consistency arises through four key choices. First is the choice of a rigid boundary condition at the outer radius of the cylinder. Second is the choice of a reduced dilatancy in the vorticity direction of the particle-stress anisotropy tensor ($\varLambda _\times < 1/2$). Together, these enable the prediction of radially inward melt segregation and compaction at outer radii, both of which are observed in experiments. The third choice is for $D_0\approx 2$, which predicts the emergence of porosity bands at 15–20$^\circ$ to the shear plane, as observed in experiments. The fourth choice is for a finite cooperativity length $\xi$, which regularises the growth-rate spectrum.

These choices are neither physically implausible nor empirically unreasonable. Therefore, we assert that laboratory experiments provide support for our hypothesis, under the conditions (i.e. the strain rate) at which they are conducted. How (indeed, if) our theory extrapolates to the much slower strain rates under natural conditions depends on the physical processes that are occurring at the grain scale. The grain-scale physics is discussed below, after we consider the relationship of our theory with that of anisotropic viscosity.

6.1. Relationship to anisotropic viscosity

A theory for anisotropic viscosity (Takei & Holtzman 2009) is also capable of explaining band angles and radial segregation (Takei & Katz 2013). The basis for this theory is a model of Coble creep, where melt provides a fast pathway for circum-grain mass diffusion. In small-strain experiments (Takei 2010), deviatoric stress causes the melt to preferentially coat grain boundaries that have normal vectors in the direction of maximum tension. The theory predicts that this anisotropy in solid-phase contiguity causes an anisotropic creep response to deviatoric stress: the deviatoric-compression direction has higher viscosity than the deviatoric-tension direction. This anisotropy gives rise to an effective dilatancy that drives radial segregation and band angle (Qi et al. 2015; Takei & Katz 2015), as also obtained here.

It may therefore make sense to think of the present, direct formulation of dilatancy as an effective description of underlying physics that is more fully described by anisotropic, Coble-creep viscosity. However, there are reasons to doubt this view. The first is that the bands that emerge in experiments on hot, partially molten rock are similar to Riedel shear zones in cold granular media (Schmocker et al. 2003; Bedford & Faulkner 2021), suggesting a common mechanism. Although both are cases of deforming granular media, Coble creep is thermally activated and does not contribute at low temperature. A second reason is that the solid-phase contiguity tensor, measured in band-producing experiments, is not suitably oriented to predict low-angle bands in viscous anisotropy theory (Takei & Katz 2013; Qi et al. 2015) (however, see Seltzer et al. 2023). And a third reason is that viscous anisotropy theory has no inherent mechanism to regularise the band growth-rate spectrum. So it may be that the granular-medium hypothesis considered here represents a distinct physical mechanism, albeit with similar implications for observable features.

Further work is required to develop experiments and analyses that can distinguish between these competing hypotheses. For example, granular physics should be tested against the hysteresis measured in oscillating stress experiments (Takei 2010). To capture this behaviour may require extension of the present theory to include a fabric tensor (Mehrabadi, Nemat-Nasser & Oda 1982) and its evolution. A more fundamental approach, however, is to develop grain-scale models that consistently integrate a set of plausible physical processes. This could clarify the conditions under which grain-boundary sliding is accommodated by dilatancy and/or Coble creep.

6.2. Physics at the grain scale

The creeping deformation of polycrystalline aggregates is known to occur by several grain-scale mechanisms: deformation of grains by the motion of lattice dislocations, shape change of grains by mass diffusion down gradients of chemical potential induced by deviatoric stress, and grain-boundary sliding whereby the centre of mass of adjacent grains moves relative to one another. This latter mechanism is typical of athermal granular flow, in which an aggregate of rigid grains is required to (locally) dilate to accommodate the geometric incompatibilities of relative motion. At higher effective stress, geometric incompatibility might instead be resolved by cataclasis or by shape change of grains through diffusion or dislocations. A sub-set of these mechanisms may simultaneously contribute to macroscopic deformation of the aggregate.

There is currently no unified model to predict the relative contributions of different mechanisms across a broad range of conditions. However, three basic expectations are relevant. First, low effective stress relative to the driving shear stress will favour dilatancy over shape change of grains. Second, higher resistance to grain–grain sliding along a grain boundary (whether viscous or associated with a frictional yield stress) will favour shape change over sliding. And third, higher homologous temperatures will favour thermally activated processes of mass diffusion and dislocation motion over cataclasis.

Saturation of the dense granular medium with a mobile, incompressible liquid may also be an important factor. (For this discussion, we consider the presence of liquid as being independent of homologous temperature even though, in the case of partial melt, the two are linked.) A greater volume fraction of liquid means that less geometric incompatibility is incurred by relative motion of grains, and hence, grain-boundary sliding is promoted. Furthermore, a greater liquid pressure relative to the mean compressive stress of the solid framework (i.e. a low effective stress) also promotes sliding. However, there are two other points to consider. First, in a poro-viscous context, non-zero effective stress causes (de)compaction, even in the absence of shear. Second, liquid transport over a finite distance through a porous medium occurs on a time scale and with a resistance controlled by the ratio of liquid viscosity to permeability. This transport is required to accommodate (de)compaction and is therefore a control on the evolution of the effective pressure.

These considerations become important in the context of dilatancy within a sealed domain of fixed volume. If the enclosed, saturated granular medium is undergoing shear, then throughout the volume there is a compressive solid stress arising from grain interactions. Dilatation can occur locally within the volume, but only if it is balanced by compaction elsewhere. There are two relevant cases. If the strain rate is nominally uniform within the domain, then dilating and compacting regions emerge by instability on length and time scales that are set internally, as in § 4.2. Alternatively, if the strain rate has an imposed gradient, then this will organise the spatial pattern of dilatation and compaction, as in § 4.1 and Appendices C and D. The particle-stress anisotropy is of fundamental importance in this latter case.

6.3. Implications for natural systems

In the shallow mantle near mid-ocean ridges, huge volumes of partially molten rock undergo deformation. These regions bear some similarity to the laboratory experiments considered here, but the strain rates are orders of magnitude smaller. Because the melt-filled pore network is vast and isolated, the effective pressure ${\rm \Delta} P$ of the solid phase may be low. Does shear cause dilatancy in this natural system? We consider this by rearranging (3.3) with ${\rm \Delta} P\approx 0$,

(6.1)\begin{equation} \frac{\mathcal{C}}{\dot{\varepsilon}_{II}} \sim \frac{D_\phi}{\zeta_\phi}, \end{equation}

where in this simple relationship, $\mathcal {C}\ge 0$ is entirely due to dilatancy. On the right-hand side is a ratio of the dilatation and compaction viscosities. Our results here suggest this ratio is ${O}(1)$ for experiments; it might be much smaller in the mantle. Hence, the question of dilatancy in a natural system may reduce to understanding how $D_\phi$ and $\zeta _\phi$ properties differ in the natural system from in the laboratory. How they scale with temperature, grain size and porosity, and whether they have a (nonlinear) dependence on strain rate are questions to be resolved by future laboratory experiments and grain-scale physical models.

At depths shallower than the mantle, crustal magmatic systems can have larger strain rates when crystal-laden melt is injected into dikes and sills (Rivalta et al. 2015). These flows have lower solid fractions, at or below jamming, and hence, more closely resemble granular suspensions (Smith 1997; Petford et al. 2020). Dilatancy should therefore be expected, and may drive crystals away from the walls of magma-filled fractures. This would reduce the effective viscosity of the magma and promote propagation.

At shallower depths and lower temperatures in the crust, seismogenic faults may experience the effects of dilatancy. The slip across faults is often accommodated by rupture of asperities or shear of a granular medium called gouge. In both cases, a compressive stress will arise (Chen & Spiers 2016). If dilatation of the pore space is possible, either by expansion of the contained air or inflow of water, there may be a decrease in friction and an unstable acceleration of slip. In contrast, if dilatation is prohibited, the compressive stress may increase friction and promote stable sliding. In either case, once sliding has ceased, viscous creep may lead to slow compaction of the gouge (or asperities) to increase contact area and harden the fault. In this case, the compaction viscosity would control the time scale for frictional state evolution (Chen, Niemeijer & Spiers 2017; Thom et al. 2023).

Indeed, the viscous constitutive law (3.3) relating compaction, dilatancy and the interphase pressure difference may be the most significant novelty of the present paper. This formulation bridges soil mechanics, suspension theory and theories for granular media with deformable grains. It may be broadly relevant in systems where deformation can occur on long time scales by irreversible creep. Such systems may be more common than is widely appreciated because slow granular processes have, until recently, gone largely unnoticed (e.g. Deshpande et al. 2021; Houssais, Maldarelli & Morris 2021).