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.


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- Reynoldsnumber 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).
However, deformation of partially molten rock inevitably includes a component of sliding along grain boundaries (e.g., Hansen et al. 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 those works, the geometric graincompatibility 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 et al. (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 • 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 • 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 et al. (2006) and Rudge & Bercovici (2015) showed that if this viscosity is (effectively) non-Newtonian with a power-law exponent of ∼6, the angle is reduced to match the observations.However, this exponent was measured by King et al. (2010) to be ∼1.5 ± 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 • (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 et al. (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 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 dilation due to shear (Reynolds 1885;Boyer et al. 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 et al. 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, ξ, of order 10× 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 et al. 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 crystalrich deforming magma (e.g., Smith 1997;Petford et al. 2020), in lower-crustal shear zones (Menegon et al. 2015), and in gouge-filled fault zones (e.g., Marone et al. 1990;Segall et al. 2010).Dilation has been considered in competition with compaction (Paterson 2001;Niemeijer & Spiers 2007) and as a microphysical mechanism responsible for rateand-state friction (Chen & Spiers 2016).It may play a role in regulating glacial sliding (Warburton et al. 2023) and in a range of geomorphological processes (Jerolmack & Daniels 2019).Dilation 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 behaviour 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.

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 ∼95% olivine grains and ∼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 ∼10 µm.Samples are hydrostatically hot-pressed to remove gas-filled bubbles prior to deformation.After hot-pressing, they have a nominally uniform melt fraction, ϕ 0 .
The experiments are conducted in a gas-medium, triaxial-deformation apparatus (Paterson 1990) with a confining pressure of 300 MPa and temperatures of ∼1225°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 1a.Cylindrical samples with height H and outer radius R are deformed by a circular platen that turns with angular velocity Ωẑ 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) φ that is axisymmetric.The velocity component U increases linearly in the ẑ and 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.
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, lowporosity lenses after a shear strain of γ ∼ 1 (Holtzman et al. 2003).These sheets are usually measured in cross-section, as in figure 1b, and hence referred to as bands.They have a characteristic spacing and form at 15-20°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°).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 center of the cylinder (Qi & Kohlstedt 2018).This migration leads to azimuthally averaged porosity, measured over the transverse section shown in figure 1c, that decreases with radius.Experiments to increasing values of total twist (reported as shear strain γ 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.

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 poroviscous theory derived by McKenzie (1984) and reviewed by Katz (2022).The melt is present with volume fraction ϕ (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 graingrain interactions during shearing flow.The constitutive law for the effective stress is then where ζ ϕ , η ϕ and D ϕ are dynamic viscosities for isotropic, deviatoric and dilational deformation, respectively, I is the identity tensor, and where are the decompaction rate, deviatoric strain-rate tensor, and particle-stress anisotropy tensor, respectively.We have introduced v s , the solid velocity field, and εII ≡ ε : ε/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 Λ 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 Λ 11 direction is taken to be the direction of flow (indicated by ∥); the Λ 22 direction is normal to the shear plane (hence we denote it Λ ⊥ ); the Λ 33 direction is the direction of the vorticity vector (and hence denoted Λ × ).The Λ 11 entry is factored out and lumped with D ϕ .Therefore, Λ ⊥ and Λ × 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 et al. 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 equation (3.1), where P j = −tr σ j /3 is the pressure of phase j.This equation states that the shearstrain rate has two possible consequences for isotropic deformation.If ζ ϕ = 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 (∆P = 0), then shear causes dilation.This has a parallel in soil mechanics, where the dilatancy angle ψ gives the kinematic relationship between shear strain and dilation (e.g., Oda & Iwashita 2020).In the present context, we can compute the dilatancy angle as tan ψ ≡ D ϕ tr (Λ) /3ζ ϕ .The more general case, of interest here, is where both the effective pressure and the compaction viscosity are nonzero.
To complete the rheological model, we require expressions for the dependency of the three viscosities on melt fraction ϕ.Empirical constraints and theoretical models of the shear viscosity η ϕ and compaction viscosity ζ ϕ 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 ϕ.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 ∼5/3 larger than the shear viscosity (Takei & Holtzman 2009;Rudge 2018).There are no empirical measurements of the dilitation viscosity D ϕ 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 ϕ to be an unknown multiple of η ϕ .Hence the viscosities are given by where η 0 is a reference value of shear viscosity at reference melt fraction ϕ 0 , λ ≈ 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 where ε∥ , ε⊥ , ε× are the normal components of the strain-rate tensor in a coordinate system aligned with simple shear.Assuming isotropic dilatancy , and into this we substitute the viscosities of equation (3.4).We find that Ψ is positive definite if 0 ⩽ D 0 < 4 5/3 ≈ 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 ξ, 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 η ϕ .We rewrite their fluidity equation in terms of a non-local viscosity, Evidently, if ξ = 0, then the non-local viscosity reduces to η ϕ .For ξ > 0, this equation imposes a minimum scale of viscosity variation.Goyon et al. (2008) measure ξ as a function of 1 − ϕ and find that it increases to about 5× the grain diameter at a solid fraction of 85%.We shall see below in §4.2 that ξ > 0 serves to regularise the spectrum of instability growth.

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 ) 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 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 section 2.

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 section 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, φ, z), as shown in figure 1a.We consider a cylinder of partially molten rock with outer radius R, azimuthal symmetry in φ and, instantaneously at t = 0, with uniform porosity ϕ 0 .At this instant, the solid flow is assumed to have zero ẑ component, a fixed azimuthal component, and an unknown radial component.This flow is described by where V (r) is the unknown radial component of the solid velocity field, Ω is the constant twist rate and H is the uniform gap between the parallel plates.We linearise the strainrate intensity under the assumption that dilatancy is driven by the forced shear such that εII ∼ Ωr/2H.(4.3) 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 εII in the radial component of force-balance equation (4.1b) to write 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 we obtain the dimensionless equation Here we have introduced 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 equation (4.6) arises from the particle-stress anisotropy tensor Λ, with the coordinates aligned such that the flow direction is φ and the vorticity direction is r.The boundary condition at the centre of the cylinder is V (0) = 0.With this constraint, equation (4.6) admits the solution 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.

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 This result demonstrates that the sign of V (r) is determined by the size of Λ × .In figure 2(a), we have chosen Λ × = 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 Λ × satisfying is also qualitatively consistent.For this range of Λ × , 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.

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, r•σ•r = −P c .Assuming that the liquid pressure is continuous at r = R, we obtain the boundary condition r 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 This condition yields a different value of A and the general solution (4.8) becomes This function is plotted in figure 2(c)-(d); the curves are computed with R = 0.3 and values of Λ × 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 Λ × 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 Λ × = 0.75, where there is a region of V < 0 at inner radii.The inner part of this region is associated with compaction C < 0, as shown by the solid curve in panel (d).The driving force is again dilatancy, but in this case with Λ × > 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.

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 ẑ, 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 where T is the torque that causes a twist-rate of Ω at t = 0. ∆F is the outward force in excess of that due to the confining pressure surrounding the sample.Using typical laboratory values for T and R in (4.14) (appendix A), we find that the excess force is on the order 10% of the force due to the typical confining pressure.
In detail, the excess force ∆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 R and Λ × .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 (eqn.(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.

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 ϕ(r) leads to a radial dependence of mobility M ϕ and aggregate viscosity η ϕ .
A series of time-dependent numerical solutions are plotted in figure 3a, coloured according to the shear strain γ at the outer radius of the domain R. Details of the numerical method are in appendix B; code is available in an online repository (Katz et al. 2023).This calculation uses R = 0.3 and Λ × = 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 C(r), shown in fig.2b.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.
In the steady state and with ξ = 0, the radial component of the solid velocity is zero and radial force-balance equation (4.1b) reduces to r • (∇ • D ϕ Λ εII ) = 0 or, after expanding and rearranging, where D ′ ϕ is the derivative of the dilation viscosity with respect to its argument, ϕ.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 Λ, is due to azimuthal dilatancy that pushes solid radially outward.The radial normal force, associated with the coefficient Λ × in Λ, has two causes: first, the gradient in radial dilatancy due to the torsional shear ( εII ∝ 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 Λ × < 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 equation (4.15) requires that D ′ ϕ < 0 independent of the specific form of D ϕ ; in other words, it requires that the dilatancy stress at fixed strain rate decreases as porosity increases.
In equation (3.4) we specified that D ϕ has the form D ϕ = D 0 η ϕ .Given the exponential dependence of η ϕ on ϕ, it follows that D ′ ϕ = −λD ϕ .With this we can solve (4.15) to give 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 Λ × ; in §5 we compare it with measurements from laboratory experiments.The logarithmic singularity in (4.16) for r → 0 is removed by non-local viscosity when cooperativity length ξ > 0. This is evident in figure 3a by comparison between the numerical solution at late time (red curve; ξ = 0.1δ) and the steady solution (black curve; ξ = 0).

Simple shear between parallel plates
Here, following the analysis of Spiegelman (2003), we investigate the stability of a twodimensional simple-shear flow with initial melt fraction ϕ 0 + ϕ 1 (x, t), where |ϕ 1 | ≪ ϕ 0 is a perturbation.A schematic diagram is shown in figure 4(a).The coordinate system is oriented such that x is in the flow direction and ŷ is in the direction perpendicular to the shear plane.We assume invariance in the ẑ (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, Chap. 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.
We use equation (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 equation (4.1b) to obtain an equation governing the solenoidal part.These are coupled to equations (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[ik(t) • x + s(t)], where 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 γ and zero compaction rate.Using the base-state solution, we solve the perturbation equations to obtain the dimensionless growth rate ṡ as a function of dimensionless wave-vector magnitude k and wavefront angle to the shear plane θ ≡ tan −1 k x /k y .In the present case, we obtain . (4.17)In this equation, ṡ has been made dimensionless by scaling with the background rate of shear γ.Wavenumber k has been made dimensionless by scaling with the compaction length δ ≡ √ 3M 0 η 0 .We have introduced the ratio ϵ ≡ ξ/δ 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 ϵ (assuming optimal orientation, θ = θ * ).The curve for ϵ = 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 ≫ 1) grow equally fast (Stevenson 1989).The use of the non-local viscosity with finite ϵ regularises the spectrum, imposing a shortwavelength cutoff at a dimensional wavelength ∼ ξ, the cooperativity scale of the non-local viscosity.Growth-rate curves in fig.4(b) have a maximum at a dimensional wavenumber k * = 1/ ξδ. (4.18) Figure 4(c) shows the growth rate as a function of the angle θ between wavefronts and the shear plane.Four curves show different values of the dilatancy viscosity prefactor D 0 ⩾ 0 (assuming Λ = 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 • , a range over which bands are subject to tension, and negative rates for angles greater than 90 • , which are subject to compression.The maximum growth rate occurs at θ * = 45 • , where band wavefronts are perpendicular to the principal tension axis.
For larger D 0 in fig.4(c), dilatancy leads to peak growth rate at low and high angles.In particular, the two maxima of growth rate occur at angles θ * 1,2 that vary with A growth-rate peak at θ * = 15 • , 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 ′ ϕ < 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 θ.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 θ.Shear localises when bands are at low or high angle to the shear plane and therefore this effect is proportional to cos 2 2θ.The combination of the two contributions of isotropic dilatancy is negative, overall, and proportional to sin 2 2θ.
Dilatancy in partially molten rock may be anisotropic, however.Experiments on dense granular suspensions, reviewed by Guazzelli & Pouliquen (2018), have obtained inconclusive estimates of Λ ⊥ , with some reporting Λ ⊥ increasing from unity with solid fraction, some reporting the opposite, and others reporting Λ ⊥ within error of unity for all solid fractions.In figure 4(d) we compare growth-rate curves as a function of band angle for Λ ⊥ = 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.

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 • , except at strains greater than about 3, at which the points spread out into a range of 10-25 • .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).
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 et al. 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., fig.4c) remains approximately independent of strain, then a forward integration of ṡ 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 Λ = I).In panel (a), for D 0 = 0, we see that without dilatancy, peak predicted amplitudes are far from observations.In panel (b), for D 0 = 2, and in panel (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 pre-factor 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 * = √ ξδ, where ξ is the cooperativity length scale and δ is the compaction length.empirically known grain diameter d and compaction length δ.The cooperativity scale ξ is understood to be proportional to grain diameter (Henann & Kamrin 2013) and hence 1/k * ∝ √ dδ.Although there is considerable uncertainty on the experimental estimates (particularly δ), a linear trend is compatible with the data.
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 normalized sample radius in figure 7. The data are colored by the magnitude of shear strain at the outer radius, which ranges from zero to ∼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.
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 2d, 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 equation (4.16).This is computed with Λ × = 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 γ ≈ 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 fig.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 ϕ 0 and λ), we could use equation (4.16) to infer Λ × .
Advancing speculatively, we assume that the experiments do approach a steady state.We note that smaller Λ × corresponds to larger dϕ/dr at steady state (fig.3b).Assuming that our estimates for ϕ 0 ≈ 0.04 and λ ≈ 27 are sufficiently accurate, the data in figure 7 are indicative of Λ × ≲ 0.4.We can estimate the timescale of porosity adjustment to steady state in the numerical solutions shown in figure 7. Referring to the porosity evolution equation (4.1c), we approximate D s ϕ/Dt by ∆ϕ/τ , where τ is the timescale over which porosity changes by ∆ϕ from its initial value ϕ 0 to its steady value at some given radius.According to (4.1c), this change is driven by decompaction at a rate (1−ϕ)C; we approximate this rate as C, which we take to be the decompaction rate at r = 0, t = 0 (c.f.fig.2b).This rate is obtained by calculating C(r) at t = 0 from the analytical solution (4.9) for radial velocity, re-dimensionalising, and evaluating at r = 0 with R ∼ 1.We then form the outer-radius strain at time τ as γ(τ ) = τ / γ = ∆ϕ/ γ C, where γ = ΩR/H is the outer-radius strain rate.This obtains where we have used (4.16) to evaluate ∆ϕ at r = 1/e 2 .For the parameters used in figure 7, this gives an outer-radius strain of γ ∼ 1 at which the simulated porosity has evolved to within a factor of 1/e of its steady value.This estimate is comparable to the numerical solution of fig. 7, but it is smaller than the empirical timescale by a factor of 10.We cannot bring these timescales into agreement by changing D 0 because this is constrained by the angle of porosity bands (fig.5), and we have already assumed that we know λ.
Reducing Λ × thus appears to be an option.According to figure 3b, reducing Λ × by a factor of 2 predicts a steady-state, radial porosity gradient much larger than observed at the largest empirically attained strain γ.Experiments to larger γ are needed to test this prediction.Alternatively, the discrepancy in timescales may be a consequence of nonlinear interaction of radial segregation with emergence of high-porosity layers, which is not captured in our models.

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 (Λ × < 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 ≈ 2, which predicts the emergence of porosity bands at 15-20°to the shear plane, as observed in experiments.The fourth choice is for a finite cooperativity length ξ, 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.

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 et al. 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.

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 shapechange 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 timescale 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.Dilation 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 dilation and compaction, as in §4.1 and appendices C and D. The particle-stress anisotropy is of fundamental importance in this latter case.

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 meltfilled pore network is vast and isolated, the effective pressure ∆P of the solid phase may be low.Does shear cause dilatancy in this natural system?We consider this by rearranging equation (3.3) with ∆P ≈ 0, where in this simple relationship, C ⩾ 0 is entirely due to dilatancy.On the right-hand side is a ratio of the dilation 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 ϕ and ζ ϕ properties differ in the natural system from in the laboratory.How they scale with temperature, grain size and porosity, and whether they have a (non-linear) 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 dilation 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 dilation 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 timescale for frictional state evolution (Chen et al. 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 et al. 2021).represents the liquid pressure acting on the plate; the third term is the indirect effect of dilatancy, which causes a net decompaction of the cylinder.
To evaluate the second and third terms in (A 6), it remains to specify V (r).We consider two cases with different boundary conditions at the outer edge of the cylinder.
A.1.No radial flow at r = R For the condition r • v s (R) = 0 and uniform porosity (e.g., at t = 0), we obtained the dimensionless solution given in equation (4.9).For this case we have V (0) = V (1) = 0 and hence the third term of (A 6) gives zero contribution.We use the solution (4.9) for dimensionless V in (A 4) to obtain This integral is evaluated by quadrature.Empirically reasonable values of the non-dimensional compaction length R are O(1).Figure 8(a) shows the non-dimensional force multiplier (in square brackets in (A 6)) for R in this neighbourhood and three values of Λ × .Depending on whether Λ × is greater than or less than 1/2, the liquid pressure (via I) contributes negatively or positively to the excess force.Recall that for this boundary-condition case, only Λ × < 1/2 gives the sign of V consistent with experiments.If Λ × ≲ 1/2, the nondimensional force multiplier is ∼ 1, meaning the the excess force is dominated by the direct effect of dilatancy.For this boundary-condition case, we therefore approximate the normal force on the parallel plates as To estimate ∆F , we take D 0 = 3 and Λ ⊥ = 1, consistent with considerations developed in the main text; we adopt values representative of laboratory experiments R = 5 × 10 −3 m and T = 10 N m (King et al. 2010).Using these, the excess outward force exerted by the sample is about 2 kN, which is about 10% of the force F c due to the experimental confining pressure (300 MPa).
A.2.No radial effective stress at r = R For the boundary condition (4.12) representing zero radial effective stress at the outer edge of the cylinder and with uniform porosity, the analytical solution for V (r) is equation (4.13).The excess axial force outward on the parallel plates is given by using this solution in equation (A 6).This case differs from that considered in §A.1 by the nonzero contribution of net decompaction, represented by the dimensionless V (1) > 0. It also differs in the contribution of the liquid pressure, associated with the integral The liquid pressure (and hence I) now depends on Λ × due to its appearance in the boundary condition.The dimensionless force multiplier (in square brackets in (A 6)) is plotted in 8(b) for three values of Λ × .Across the range of R and Λ × considered, the force multiplier ranges from about 20% to about 60%.Dilatancy in the radial direction evidently reduces the outward normal force on the parallel plates.Using the experimental and theoretical values quoted in §A.1, the expected outward excess force is about 1 kN, which is about 5% of the force due to confining pressure.

Appendix B. Finite-time and steady models of parallel-plate torsion
By considering the system of conservation equations (4.1), we can derive a model for the evolution of the radial distribution of porosity in parallel-plate torsion.We again assume the simplified, axisymmetric torsional flow considered in §4.1, This flow is substituted into the compaction equation (4.1a), which is then integrated subject to boundary conditions V (R) = 0 and ∂P ℓ /∂r| R = 0.The radial component of the momentum conservation equation (4.1b) is simplified with (B 1) and used to eliminate the pressure gradient.We non-dimensionalise the result, along with equation (4.1c) and equation (3.6) using characteristic scales where the dimensionless local viscosity is η ϕ = exp[−λ(ϕ − ϕ 0 )] and ϵ ≡ ξ/δ.It is important to note that the dilation-viscosity coefficient D 0 does not appear in these equations; it appears only in the relationship between strain at the outer radius γ and dimensionless time t, t = γD 0 /6.(B 4) This relationship indicates that a given dimensionless time (and hence amount of porosity change) is reached at smaller outer-radius strain when D 0 is larger.In other words, increasing D 0 promotes radial melt segregation.
Boundary and initial conditions imposed on the system (B 3) are This system (B 3) with conditions (B 5) can be solved numerically by one-dimensional finite-volume discretisation on a grid with uniform intervals in r and t.Velocities are stored at nodes (the points that connect intervals); porosity and viscosity are stored at interval centres.The code, developed in the framework of the Portable, Extensible Toolkit for Scientific Computation (PETSc, Balay et al. 2023Balay et al. , 1997)), is available in an online repository (Katz et al. 2023).
Numerical solutions indicate that ϕ(r, t) tends toward a steady state in which V (r) = 0.This state can be determined analytically for ϵ = 0 and hence when η ϕ = η ϕ .Then, (B 3a) reduces to 1 Solving and rewriting in terms of the porosity gives where we have determined the constant of integration using global conservation of liquid mass, 2 1 0 rϕ dr = ϕ 0 .

Appendix C. Cone-and-plate torsional flow
Cone-and-plate torsional flow is another geometry that has been used in experiments on dense granular suspensions (Guazzelli & Pouliquen 2018).Although there are currently no deformation experiments on partially molten rock in this configuration, we consider the predicted radial flow for completeness.We use the flow described by (4.2) but now take the gap H to be increasing linearly with radius, H = r.With this choice, the shearstrain rate is independent of radius.As a consequence, the non-dimensional governing equation becomes ∂ ∂r where we have scaled the radial velocity component V with In the large-compaction-length limit of R ≫ 1, the radial velocity solution is asymptotic to V (r) = (Λ × −1)(r ln r)/2.A small-R matched asymptotic solution exists but converges only for very small R.
Figure 9(a) shows plots of solution (C 4) for three values of R and for R → ∞.We have chosen Λ × = 0.45 for plotting purposes.The pattern of radial flow is similar to the parallel-plate case, but with a maximum speed that is larger and shifted toward the centre of the cylinder.The associated decompaction rate (panel (b)) increases sharply near r = 0.This decompaction is balanced by weak compaction near r = 1.
For solutions as in fig.9(a) with no flow at the outer boundary, we require Λ × < 1 to achieve a dimensional V > 0 consistent with experiments.This is less restrictive than the constraint (4.10) obtained from parallel-plate torsion with the same boundary conditions.
The zero-effective-stress boundary condition r • σ eff • r = 0 at the outer boundary, after non-dimensionalising with [V ] from eqn. (C 2), is written The analytical solution in this case is (C 6) Plots of this solution with R = 0.3 are shown in figure 9(c).The radial pattern of flow and compaction are different from the zero-solid-flow boundary condition.It is important to note that for the zero-stress boundary condition, the compaction rate need not integrate to zero over the domain because the flow (solid and liquid) at the outer boundary is non-zero.
The results in this appendix demonstrate that the radial profile of flow is sensitive to the geometry of deformation, as anticipated from previous work (e.g., Morris & Boulay 1999).They also show that the radial flow is sensitive to the outer boundary condition.These results are valid at t = 0, when the porosity, permeability, shear viscosity, and dilatancy viscosity are uniform.At later times, melt segregation causes spatial variations in porosity and hence in these coefficients.

Appendix D. Poiseuille flow through a pipe
Here we consider Poiseuille-like flow of partially molten rock along an infinite, straight pipe with circular cross-section and radius R. At t = 0, the porosity is uniformly ϕ 0 .A cylindrical coordinate system (r, φ, z) is aligned with the the axis of the pipe.The radial direction is shear-plane perpendicular and the azimuthal direction aligns with the vorticity.An axisymmetric flow is driven by a pressure gradient −G in the z direction.We assume translational invariance in z.The solid velocity, compaction rate and deviatoric strain-rate tensor are then Here, again, R ≡ √ 3η 0 M 0 /R is the ratio of the compaction length to the outer radius.For a cylindrical domain with azimuthal symmetry that extends inward to r = 0, the radial velocity must vanish on the axis, V (0) = 0.The rigid outer wall requires that

Figure 1 .
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 (ϕ0 = 0.04, γ = 1.5;King et al. 2010).(c) A transverse section showing radially inward migration of melt (ϕ0 = 0.10, γ = 5.0; Qi et al. 2015).Cracks visible in panels (b) and (c) are a consequence of the rapid quench and decompression after deformation.
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 η ϕ in equation (4.1b), which couples it to equation (3.6) governing the non-local viscosity.For consistency, we use the non-local viscosity in (3.4) to compute non-local compaction ζ ϕ and dilitation D ϕ viscosities.The second modification is the last term on the right-hand side of equation (4.1b), which captures the hypothesised dilatancy effects.The classical model is recovered for ξ, D 0 → 0.

Figure 4 .
Figure 4. Growth rate of sinusoidal perturbations under a simple-shear flow from eqn. (4.17).(a) Schematic diagram showing a finite region of the infinite domain.Grayscale shows the perturbed porosity field.(b) Growth rate as a function of wavenumber k for D0 = 0 at θ = 45 • .Circles represent the growth rate computed at k * = ϵ −1/2 .(c) Growth rate angular factor as a function of θ with Λ ⊥ = 1, and values of D0 given in the legend.(d) Growth rate angular factor as a function of θ with D0 = 2, and values of Λ ⊥ given in legend.

Figure 5 .
Figure 5. Angle spectra of porosity-band amplitude as a function of shear strain γ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 ṡ(t) from equation (4.17) with Λ ⊥ = 1 and dimensionless k(t = 0) = k * = ϵ −1/2 .Each panel has a different magnitude of dilatancy: (a) D0 = 0; (b) D0 = 2; (c) D0 = 3.

Figure 6 Figure 6 .
Figure6.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).

Figure 7 .
Figure 7. Radial distribution of porosity ϕ(r) normalised by the initial porosity ϕ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 (γ = 0) experiment.Colours represent the shear strain at the outer radius γ(R).The black, dotted line represents the steady-state solution (4.15) with Λ× = 0.4, ϕ0 = 0.04, λ = 27.Dashed curves are numerical solutions to the system (B 3) with the same parameters as for the steady curve and also D0 = 2, R = 0.3.The numerical solution is plotted at finite values of γ(R), as given by their colour.Appendix B gives details of the numerical method.

[
with inner boundary condition V (0) = 0 has the general solutionV (r) = R(Λ × − 1) [K 1 (r/R) − R/r] + BI 1 (r/R), (C 3)where I n (z) is the modified Bessel function of the first kind and K n (z) is the modified Bessel function of the second kind.The solution with outer boundary condition V (1) = 0 is
v s = V (r)r + W (r)ẑ, subject to a symmetry condition at r = 0 and a no-slip condition at r = R gives the standard Poiseuille solution W (r) = G R 2 − r 2 /4η 0 .As we did for torsion, we linearise by neglecting the contribution of dilatant flow to the second invariant of the strain rate and obtain εII ∼ 1 2 |∂W/∂r| = Gr/4η 0 .We use this in the radial component of force-balance equation (4.1b) to write 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 (2Λ ⊥ − Λ × ).(D 5)