Consequences of viscous anisotropy in a deforming, two-phase aggregate. Why is porosity-band angle lowered by viscous anisotropy?

Consequences of viscous anisotropy in a deforming, two-phase aggregate. Why is porosity-band angle lowered by viscous anisotropy?


In laboratory experiments that impose shear deformation on partially molten aggregates of initially uniform porosity, melt segregates into high-porosity sheets (bands in cross-section). The bands emerge at to the shear plane. A model of viscous anisotropy can explain these low angles whereas previous, simpler models have failed to do so. The anisotropic model is complex, however, and the reason that it produces low-angle bands has not been understood. Here we show that there are two mechanisms: (i) suppression of the well-known tensile instability, and (ii) creation of a new, shear-driven instability. We elucidate these mechanisms using linearised stability analysis in a coordinate system that is aligned with the perturbations. We consider the general case of anisotropy that varies dynamically with deviatoric stress, but approach it by first considering uniform anisotropy that is imposed a priori and showing the difference between static and dynamic cases. We extend the model of viscous anisotropy to include a strengthening in the direction of maximum compressive stress. Our results support the hypothesis that viscous anisotropy is the cause of low band-angles in experiments.

1 Introduction

In laboratory experiments, forced shear deformation of nominally uniform, partially molten rocks causes melt segregation into high-porosity bands oriented at low angle () to the shear plane (Holtzman et al., 2003; Holtzman and Kohlstedt, 2007; King et al., 2010). Stevenson (1989) predicted the emergence of such bands in a self-reinforcing feedback arising from the porosity-weakening of the crystalmagma aggregate, but the angle predicted by this theory was (Spiegelman, 2003), much higher than observed. The low angle of high-porosity bands is widely thought to provide an additional constraint on the rheology of the aggregate, but it has proven challenging to understand. Katz et al. (2006) found that non-Newtonian viscosity with a high sensitivity to stress could reproduce the low angle of bands, but King et al. (2010) subsequently showed that the viscosity of experiments that produce low-angle bands is actually close to Newtonian.

Theory by Takei and Holtzman (2009a, b) of anisotropic viscosity under diffusion creep of a partially molten aggregate represents a possible solution. This theory is motivated by observations of the coherent alignment of melt-pockets between solid grains under a deviatoric stress (e.g. Daines and Kohlstedt, 1997; Zimmerman et al., 1999) and of the enhancement of diffusion creep by melt at grain boundaries and triple junctions (e.g. Cooper et al., 1989). The melt is a fast pathways for diffusional transport of solid constituents around grains; the alignment of melt with respect to the principal-stress directions hypothetically results in anisotropic viscosity of the aggregate (Takei and Holtzman, 2009a).

Analysis of the theory of anisotropic viscosity by Takei and Holtzman (2009b), Takei and Katz (2013), Katz and Takei (2013), and Allwright and Katz (2014) shows that it introduces qualitatively different behaviour from previous models with isotropic (and even power-law) viscosity. Shear and normal components of stress and strain-rate are coupled under viscous anisotropy; as a result of this coupling, a gradient in shear stress becomes a driving force for melt segregation that is not present in the isotropic system. Under Poiseuille flow, melt segregates toward higher-stress regions; under torsional flow, compressive hoop stresses drive the solid outward and the magma inward. The mechanics of this “base-state” melt segregation are explained in detail by Takei and Katz (2013). An experimental test of radial melt segregation in torsional flow by Qi et al. (2015) shows striking consistency with predictions.

Furthermore, theoretical work has demonstrated that there is a connection between the strength of anisotropy and the angle of high-porosity bands that emerge by unstable growth. This was shown with linearised stability analysis (Takei and Holtzman, 2009b; Takei and Katz, 2013) and numerical simulations (Butler, 2012; Katz and Takei, 2013) where the strength and orientation of anisotropy are assumed to be known and are imposed a priori. In those static-anisotropy calculations, high-porosity bands emerge at low angles to the shear plane only when viscous anisotropy is at or near saturation. This is a rather restrictive condition that may be incompatible with the robust appearance and consistently low angle of bands in experiments (Holtzman and Kohlstedt, 2007). However, in numerical simulations that allow anisotropy strength and direction to vary dynamically in space and time (Katz and Takei, 2013), band angles are significantly lowered and appear to be less sensitive to the mean strength of anisotropy. These findings raise several basic, unanswered questions: Why do the mechanics of viscous anisotropy give rise to low-angle bands? Why is dynamic anisotropy more effective in this regard than static anisotropy? What are the general conditions under which low-angle, high-porosity bands should form?

The present manuscript addresses these questions through a combination of linearised stability analysis and physical reasoning. The crucial, enabling advance is to perform the analysis in a coordinate system that is rotated to align with the porosity bands (rather than with the plane of shear). This drastically simplifies the expressions for growth rate under static anisotropy (Takei and Katz, 2013), making them readily interpretable in physical terms. Moreover, it allows us to extend the analysis to dynamic anisotropy in a form that exposes the physical differences from static anisotropy. Finally, the same coordinate rotation clarifies the physical reason for low angles under isotropic, non-Newtonian viscosity.

The manuscript is organised as follows. In the next section, we briefly discuss the nondimensionalised governing equations and present an anisotropic, viscous constitutive model for the two-phase, partially molten aggregate. The full, non-linear system is solved numerically in §3 for static and dynamic cases, to elucidate the questions listed above. The coordinate rotation is introduced and the linearised stability analysis is developed in §4. In particular, §4.3 develops an expression for the growth-rate of porosity perturbations under the fully dynamic model of §2. This expression is challenging to understand and so we subsequently consider it under reducing assumptions of static anisotropy (§5.1), which includes the simplest case of Newtonian, isotropic model. We build on this to explain the full complexity in §5.2 and §5.3. We conclude with a summary and discussion of the results in terms of the motivating questions.

2 Governing and constitutive equations

In the theory of magma/mantle interaction, the macroscopic behaviour of a two-phase aggregate is treated within the framework of continuum mechanics (e.g. Drew, 1983; McKenzie, 1984). This theory is concerned with the evolution of macroscopic fields including the volume fraction of melt or porosity , the velocity of the solid phase , the liquid pressure (compression positive), and the bulk or phase-averaged stress tensor where is the stress tensor of the solid phase (tension positive). Further details of the two-phase-flow theory were previously presented (e.g. Takei and Katz, 2013; Rudge et al., 2011) and are not repeated here.

We proceed directly to the nondimensional governing equations,


and refer the reader to Takei and Katz (2013) and references therein for details of the derivation and rescaling. In the system (1) we have introduced the differential stress tensor . Also, we have excluded body forces and assumed that the permeability of the solid matrix is a function of the porosity only, proportional to , where is a reference porosity and is a constant. is the nondimensional compaction length and is a rheological parameter explained below. To close the system, a constitutive relationship that relates the differential stress and strain rate is required.

Takei and Holtzman (2009b) and Takei and Katz (2013) proposed a model of anisotropic viscosity caused by stress-induced microstructural anisotropy. In partially molten rocks, the melt phase is contained within a permeable network of tubules between grains. The solid matrix is formed by a contiguous skeleton of solid grains. The area of grain-to-grain contact is known as the contiguity. Contiguity is the microstructural variable that determines the macroscopic (i.e., continuum) mechanical properties of the matrix (Takei, 1998; Takei and Holtzman, 2009a). Although the equilibrium microstructure developed under hydrostatic stress has isotropic contiguity, deviations from the equilibrium microstructure have been observed in experimentally deformed, partially molten samples (e.g. Daines and Kohlstedt, 1997; Takei, 2010). Based on these observations, we infer that under a differential stress, the grain-to-grain contacts with normals that are parallel to the maximum tensile stress () are reduced in area; similarly, the areas of those with normals parallel to the maximum compressive stress () are increased. Using a microstructure-based model of aggregate viscosity (Takei and Holtzman, 2009a) and a coordinate transformation (Takei and Katz, 2013), the constitutive law and the viscosity tensor are


For simplicity, we consider a two dimensional problem, in which the plane is parallel to the plane. Therefore only the two-dimensional version of is written in (2b). Only 6 of the 16 components are shown due to the symmetry of under the exchange of and , and , and and .

The factor in front of the matrix represents the normalised shear viscosity ; it decreases exponentially with increasing melt fraction , and so is called the porosity-weakening factor. We take based on the experimental results (e.g. Mei et al., 2002). The parameter represents the bulk-to-shear viscosity ratio, , which is assumed to be constant (=5/3) based on theoretical results by Takei and Holtzman (2009a) (although see Simpson et al., 2010a, b). Parameters , , and represent the magnitude and direction of microstructural anisotropy: and quantify the amplitude of contiguity reduction and increase, respectively, in the directions of principal stress and ; represents the angle that the most-tensile stress () direction makes with the -axis of the coordinate system. Using the local differential stress , is given by


and and are modelled as


where represents the amplitude of deviatoric stress. The detailed forms of the functions in (4) are poorly constrained, owing to a lack of experimental data. The form of was chosen based on the constraints that is less than or equal to 2 (Takei and Katz, 2013) and increases with increasing differential stress (Daines and Kohlstedt, 1997; Takei, 2010). The parameter is assumed to be a constant that is probably between 0 and 1. In the present study, is parameterised by and , which control the stress-offset and slope of increase. In Takei and Katz (2013), and was parameterised by alone. Parameters and are newly introduced here.

For simplicity in previous linearised analyses (Takei and Katz, 2013; Allwright and Katz, 2014), parameters , , and were fixed to their initial values. We call this simplifying assumption the static anisotropy model. In contrast, the complete model with stress-dependent direction and magnitude is called the dynamic anisotropy model. Katz and Takei (2013) discovered a remarkable difference between static and dynamic anisotropy; this difference motivates the present study and is demonstrated in the next section.

3 Numerical solutions

Numerical solutions of equations (1) and (2) highlight the difference between results obtained for static and dynamic anisotropy. The solutions are computed with a finite-volume method on a fully staggered grid that is periodic in the -direction; in this section, the axis is taken parallel to the initial flow direction. No-slip, impermeable boundary conditions enforce a constant displacement rate of plus or minus on the top and bottom boundaries, respectively. A semi-implicit, Crank-Nicolson scheme is used to discretise time and the hyperbolic equation for porosity evolution is solved separately from the elliptic system in a Picard loop with two iterations at each time-step. The solutions are obtained in the context of the Portable, Extensible Toolkit for Scientific Computation (PETSc, Balay et al., 2001, 2004; Katz et al., 2007). Full details and references are provided by Katz and Takei (2013).

Figure 1: Comparison of numerical solutions to (1) and (2) with static and dynamic anisotropy. For both calculations, and the domain is discretised into 600300 square cells. (a) Porosity field at a strain of 0.75 for a simulation with and throughout the domain. (b) Porosity field at a strain of 1.25 for a simulation with anisotropy calculated according to (3) and (4) with . (c) Spectral power binned by wavefront-angle to the shear plane (after Katz et al., 2006) for the porosity fields shown in panels (a) and (b). Each spectrum is normalised by its maximum power. (d) and (e) Two-dimensional histograms derived from the simulation with dynamic anisotropy at a strain of 1.25 (after Katz and Takei, 2013, Fig. 12). Red dashed lines have a slope given by the ratio of perturbation quantities and from the stability analysis in §4.

Figure 1 compares solutions with fixed and dynamic anisotropy. In panel (a), anisotropy parameters are prescribed as , ; in panel (b), these parameters are calculated cell-wise using equations (3) and (4), with . Both calculations have (compaction length equal to domain height) and are initialised with the same porosity field, , where and . is a smooth, random field with unit amplitude, generated by filtering grid-scale white noise to remove variation at wavelengths below 15 grid-cells. Because the growth-rate of porosity perturbations differs for fixed and dynamic anisotropy, the simulations are shown at different values of the average simple-shear strain .

The different orientation of high-porosity features is evident in panels (a) and (b): dynamic anisotropy is associated with lower angles. This is quantified by the power spectrum in panel (c), where the power from a 2D fast-Fourier transform of the porosity field is binned according to the angle between the wavefront and the shear plane (Katz et al., 2006). Dynamic anisotropy produces a peak at whereas static anisotropy produces a peak at . There is also a high-angle () peak for static anisotropy (corresponding to features visible in panel (a)) that does not survive at large strain. Panels (d) and (e) show the covariation of and with in panel (b); black dotted lines indicate mean values. These means are closely matched with the parameter values used in the fixed-anisotropy simulation. It is therefore clear that the difference in the dominant band angle (panel (c)) arises from the coupling between stress and the variations in and . What is unclear, however, is the physical explanation for this difference and, indeed, why viscous anisotropy gives rise to bands at angles less than to the shear plane at all. We clarify these points below.

4 Linearised analysis with a perturbation-oriented coordinate system

Let be the initial amplitude of a porosity perturbation. We express the problem variables as a series expansion about the base-state in which the porosity is uniform and equal to . We truncate the series after the first-order terms,


The first term of (5) with index represents a simple-shear flow and its associated anisotropy, which is the base-state solution of order one (), corresponding to the uniform porosity . The second term of (5) with index 1 represents the perturbation of order caused by . By substituting (5) into equations (1), using , and balancing terms at the order of , we derive the governing equations for the perturbations as


where .

Figure 2: Schematic diagrams of the coordinate axes and porosity perturbation. (a) The coordinate system for the linearised analysis is taken such that the axis is parallel to the initial perturbation wavefronts. The shear plane of the base-state, simple-shear flow is then rotated by an angle . (b) The base-state normal stress is oriented parallel to the initial perturbation wave vector. (c) The base-state shear stress is parallel to the wavefronts.

Following previous studies, the porosity perturbations take the form of a plane wave oriented at a given angle to the base-state shear plane. Past workers chose to align the coordinate system with the base-state shear plane, such that the base-state strain rate tensor has a simple form (e.g. Spiegelman, 2003). Although the coordinate system was so aligned in the numerical models above, in this section the coordinates are rotated such that the axis is parallel to the wave-vector of the initial perturbation, as shown in Figure 2. With this choice, again represents the angle between the perturbation wavefronts and the base-state shear plane. However, in the rotated coordinate system, depends on both the direction of and the orientation of the bands. The base-state direction of maximum tensile stress makes an angle to the shear plane (eqn. (3)) and so for a coordinate rotation by , we have (Figure 2a).

In the following part of this section, we give an outline of the linearised approach, which shows that the new coordinate system reduces the complexity of the analysis and exposes the physical mechanisms of perturbation growth. This enables us to clarify the mechanics leading to low-angle bands in complicated problems such as under dynamic anisotropy. We first consider the base-state, simple shear flow at the order of 4.1) and then the linearised governing equations at the order of 4.2). Finally, in §4.3, we obtain the growth rate of porosity perturbations for the most general case of dynamic anisotropy. The result obtained is used in §5 to clarify the mechanisms of low-angle-band formation.

4.1 Base-state simple shear flow

Using the angle between the initial perturbation wavefronts (aligned with the direction) and the base-state shear plane (Fig. 2a), the components of the base-state strain rate tensor in the rotated coordinate system are


As shown in Figs. 2b and 2c, and represent, respectively, the base-state tensile and shear stresses normal and parallel to the perturbation wavefronts, which play important roles in understanding the growth of these perturbations. Noting that , these components are given by


with . In Figure 3, and are plotted as a function of angle . An understanding of their systematics is needed to interpret the results of the stability analysis.

Figure 3: Normal stress (tension positive) and shear stress produced by the base-state simple shear flow as functions of angle between the shear plane and the -axis (Fig. 2a). We use here and throughout (Takei and Katz, 2013). Each panel is computed with a different set as labelled above and to the right of the panels.

Considering first the solid curves representing normal stress, we see that for (isotropic, top-left panel), follows the expected pattern of . It is tensile for and compressive for (Spiegelman, 2003). However, as increases (for , top row of panels), becomes negative (compressive) at all angles. The mechanism for this change is two-fold. First, decreases the viscous resistance to extension in the direction and reduces the maximum tensile stress. This is because increasing reduces for angles near (Figure 4a, black solid line). Superimposed on this is a compressive stress around and that emerges as a consequence of shear strain rate coupled to normal stress via the viscosity (Figure 4a, gray solid line). The product is negative for all , motivating us to name this coupling “shear-strain-induced compression.” This non-trivial result comes from the fact that the stress-induced softening occurs in the tensile () direction, as schematically illustrated in Takei and Katz (2013) Figure 7b.

The effect of increasing on is shown down columns in Figure 3. Similar to , couples the shear strain rate to compressive normal stress for angles near and via the viscosity. In contrast to , however, strengthens the aggregate in the direction and increases the normal stress amplitude near (Figures 3 and 4b). As a result, the sign-change of normal stress caused by occurs in a limited range of and (Fig. 3, bottom-left panel).

The shear stress curves in Figure 3 (dashed lines) also change with increasing and/or . For zero anisotropy in panel (a), the shear stress follows , as expected for coordinate rotation only. Anisotropy does not change the mean of , as required by symmetry of the stress tensor. Increasing (or ) has an overall weakening (or strengthening) effect, changing only the amplitude of . This is in contrast to the effect of anisotropy on normal stress, which has a strong dependence on angle .

Figure 4: Viscosity components from equation (8) as functions of angle between the shear plane and the -axis. The curves are computed for the anisotropy parameters (a) and , and (b) and . The thin black lines represent each component for the isotropic case ().

4.2 Growth of porosity perturbations

For simplicity in this linearised analysis, we consider the case of liquid viscosity , giving a non-dimensional compaction length . In this limit, liquid segregation over any length scale occurs at vanishingly small pressure gradients. Therefore, the pressure gradient terms in equations (6c) are negligible. Pressure then drops out of the problem and we no longer need to solve equation (6b) (which has become indeterminate anyway!). This is equivalent to considering only the subset of perturbations with wavelengths much smaller than the dimensional compaction length (e.g. Katz et al., 2006).

Using the initial wavenumber vector , the porosity perturbation at time is


which accounts for rotation of the wave-vector due to advection by the base-state flow (Spiegelman, 2003). At , by choice of the coordinate system, perturbations are uniform in the direction. Therefore, partial derivatives of the first-order quantities with respect to are zero. Using , , and , equations (6) become


Equations (10b) and (10c) are obtained after an integration in the direction; boundary conditions are not needed because the domain is infinite and the first-order fields are periodic.

The right-hand side of equations (10b) and (10c) represent and , respectively. Since pressure gradients are negligible, these stresses must be spatially uniform for the system to be in balance. Therefore the first-order product of viscosity and strain rate must sum to zero; viscosity reduction associated with porosity perturbations (within the first square brackets in the right hand side) is compensated by the strain rate perturbations (within the second square brackets).

To facilitate the physical interpretation of equations (10b) and (10c), these equations are re-expressed as


with equivalent (“forcing”) stresses


Equation (11) relates the strain-rate response of the system to the forcing stresses defined by (12). From equation (11), normal strain rate in the direction (the component that is most relevant to the perturbation growth) can be expressed as


where and are the compliances defined by


The forcing stresses, and , are not externally applied (like those causing simple shear), nor are they the first-order stress perturbations and (these are both equal to zero). Instead, they are equivalent stresses that are created internally as a consequence of the base-state flow acting on the viscosity change associated with porosity perturbations. Moreover, under dynamic anisotropy, these forcing terms also depend on the strain-rate perturbations and hence equation (13) does not always give an explicit solution for . Nonetheless, equation (13) enables us to separate the mechanics into two simpler parts: the forcing, and , and the compliance, and , where the latter represents the system response to forcing with unit amplitude. This decomposition is helpful to understand the detailed (and rather complicated) mechanisms of the different models considered here.

4.3 General solution

Equations (10) are solved here to obtain an explicit expression for for the full model of dynamic anisotropy. The first-order viscosity tensor is written in terms of the porosity and anisotropy perturbations , , , and . The anisotropy perturbations are then expressed in terms of the porosity and strain-rate perturbations , , and . These calculations are sketched in Appendix B. The components of the first-order viscosity tensor are then substituted into equations (10b) and (10c), which are manipulated to solve for and as functions of . The normal strain rate obtained by this approach is substituted into (10a) to give an expression for the growth-rate of perturbations,


with the compliances defined by (14) and dynamic factors




The constant coefficients and that express the sensitivity of the growth rate to dynamic anisotropy perturbations are


where and are defined by equations (34) and (36), respectively. We do not attempt to physically interpret the detailed form of and . It is important to note, however, that for the static anisotropy model, and are zero, and .

5 Physical interpretation in various limits

The growth rate in eqn. (15) is a general result for the full model presented in §2 above (with the sole assumption of ). To build up a physical understanding of this equation, we return to the simpler case of static anisotropy, which includes the simplest case of the Newtonian, isotropic model. The static anisotropy model has previously been studied by Takei and Katz (2013) using linearised analysis. However, the mathematical complexity of their results precluded a detailed mechanical interpretation. A reconsideration using the perturbation-oriented coordinate system enables a physical understanding of the instability mechanism and the rheological control on the dominant band angle. These are needed to understand the more complicated, dynamic model. To facilitate this (in §5.15.2), we make the simplifying assumption that  — that there is no contiguity increase in the -direction. The effect of non-zero is discussed in section 5.3, where we show that its role is minor compared to that of , the contiguity decrease in the -direction.

5.1 Static anisotropy

When , the mechanical equilibrium conditions (10b) and (10c) are written as


Comparison with equation (11) shows that the forcing stresses are given by and . These forcing stresses are caused by the base-state tensile and shear stresses acting on the porosity perturbation by way of porosity weakening rheology () as depicted in Figs. 2b,c. In this simple model, and are given in terms of the porosity perturbation , and hence eqn. (13) provides an explicit solution for .

It is evident from (10a) that the normal strain rate causes an increase in the amplitude of porosity perturbations; the shear strain rate does not cause the porosity to change. When the viscosity is anisotropic, both the normal and the shear stress drive and hence contribute to perturbation growth. This does not occur under isotropic viscosity.

Instability mechanism in the isotropic system

For an isotropic aggregate (), in equations (11) and (19), and the compliance that couples shear stress to normal strain rate is zero. In this case, is driven only by the base-state normal stress . The growth rate is given by


Perturbations are unstable under tensile stress () normal to the perturbation wavefronts, stable under compressive normal stress (), and unaffected by shear stress . We therefore term this the tensile-stress-induced instability, or tensile instability. When band angle relative to the simple shear flow is , the tensile stress attains its maximum (Figure 3a) and hence the growth rate is also at a maximum, as shown in Figure 5a and by Spiegelman (2003). The occurrence of the tensile instability in a porosity-weakening, two-phase aggregate was first predicted by Stevenson (1989).

Two instability mechanisms in the anisotropic system

For an anisotropic aggregate ( and/or ), there is a coupling between shear and normal components via . In this case, is forced by both normal stress across perturbations and shear stress along perturbations. The growth rate is


using the compliances given by equations (14). The first term on the right-hand side of (21) represents the tensile instability, generalised to the anisotropic aggregate. The 2nd term represents a shear-stress-induced instability that does not occur in the isotropic system. The total growth rate versus band angle is plotted in the top row of panels of Figure 5 for various anisotropy amplitudes (thick lines). Consistent with previous work, as is increased, the single growth rate peak splits into two peaks at low and high angles to the shear plane (Fig. 5c). Because the lower-angle peak dominates the higher angle peak after a finite time (Katz et al., 2006; Takei and Katz, 2013), this result means a significant lowering of the dominant band angle by the viscous anisotropy — if the magnitude of anisotropy is sufficiently close to saturation ().

Figure 5: Characteristics of the static anisotropy model as a function of the angle between porosity perturbations and the shear plane. Each column is for a different value of . In all cases, , , and . Top row. Growth rate of perturbations from eqn. (21). The heavy line represents the total growth rate; light lines show the growth rate decomposed into two parts: the contribution from the tensile instability (1st term of eqn. (21)) and that from the shear instability (2nd term). Bottom row. Compliances and in eqn. (21).

How viscous anisotropy causes lowering of band angle

Although the effect of viscous anisotropy is evident from the total growth rate shown in the top panels of Figure 5, it is not immediately obvious why the dominant band angle is lowered by viscous anisotropy. The physical mechanism can be understood by considering the tensile and shear components of the growth rate independently (1st and 2nd terms of (21), respectively). In the top row of Figure 5, these two growth rates are plotted separately for various values of (thin solid curve for tensile instability; thin dashed curve for shear instability). Comparison of panels (a), (b), and (c) reveals that the peak split occurs through (i) stabilisation of the tensile instability and (ii) emergence of the shear instability with increasing magnitude of anisotropy . We consider each of these in turn.

To understand why viscous anisotropy stabilises the tensile instability, we return to the systematics of the base-state stress (section 4.1). Comparison of the three columns of Fig. 3 shows that as increases, the tensile stress decreases in amplitude and becomes compressive at all angles. With , the first term in equation (21) is always less than or equal to zero, and hence stable.

To understand why viscous anisotropy destabilises the shear mechanism, we consider the coupling between the shear stress that drives the instability and the normal strain-rate that is responsible for its growth. As shown by equation (13) with and , the shear stress is coupled to the normal strain rate via . The angular dependence of is shown by dashed curves in the bottom-row panels of Figure 5. If is positive, then is positive (or is in phase with ) and the shear mechanism contributes to unstable growth of porosity perturbations. In fact, this product is positive (or zero) for all , enabling us to name this coupling “shear stress-induced expansion.” This non-trivial result comes from the assumed microstructural behaviour: that stress-induced softening occurs in the tensile () direction, as illustrated in Takei and Katz (2013) Figure 7a. So the porosity perturbation grows because of the shear mechanism, for which the low angle is favorable.

Summary of static anisotropy model

As a recap and summary, note that under isotropic viscosity, the growth of bands at to the shear plane is caused by a tensile instability (Stevenson, 1989; Spiegelman, 2003). In contrast, under anisotropic viscosity, the peak growth rate of bands is controlled by a distinct shear instability. Although the peak growth rate of the shear instability occurs at , stabilisation at these low angles by the tensile mechanism acts to give a maximum in the combined growth rate at .

The comparison between isotropic and anisotropic systems developed above is summarised in the first two rows of Table 1. The tensile instability is separated into porosity-weakening , which is fundamental to all models, and the tensile stress across bands , which affects both isotropic and anisotropic cases. A shift of to more negative, compressive values (represented by ) stabilises the tensile instability. In contrast, the difference in shear instability can be simply shown by its existence or non-existence ( or –). It is the leading-order terms and that are responsible for these differences.

Katz et al. (2006) extended the analysis of isotropic viscosity to include a power-law dependence of viscosity on strain rate (or, equivalently, on stress). They showed that strain-rate weakening viscosity leads to lowering of band angle. In the discussion section, we compare the angle-lowering mechanism of viscous anisotropy to that of the power-law viscosity. This is enabled by a reanalysis of the power-law model using the rotated coordinate system.

(tensile) (shear)
Model additional factor Dominant angle
Isotropic Newtonian
Static anisotropy
Dynamic anisotropy
Isotropic power-law
Table 1: Summary of band formation models. : exists; –: does not exist; : is modified.

5.2 Dynamic anisotropy

We return to the full expression for the growth rate of bands, eqn. (15), to develop a physical understanding of why dynamic anisotropy lowers band angles, as observed in the numerical solutions (Fig. 1). To do so we take and again make the simplifying assumption that (though see §5.3, below).

The perturbations in and under dynamic anisotropy are obtained by linearisation of equations (3) and (4) with respect to the stress perturbation . The expansion is conducted around the base-state values . In Appendix B, we show that the sensitivity of to variations in deviatoric stress is given by the parameter


This parameter allows us to write , and hence to see that static anisotropy corresponds to the case where . The situation for is more complicated because there is no single parameter that controls its sensitivity to deviatoric stress; variations of can either be fully considered or fully neglected. Fortunately, numerical and analytical results show that these variations () play an insignificant role in the understanding of band angles, and hence we need consider only the magnitude of anisotropy . This is achieved by looking at the dependence of key quantities (especially ) on .

Figure 6: Characteristics of the dynamic anisotropy model for various values of as a function of the angle between porosity perturbations and the shear plane. In all panels, , , , and . (a) Full growth-rate from eqn. (15). (b) Growth-rate from eqn. (15) split into the tensile-instability term (solid lines) and the shear-instability term (dashed lines). (c) Band-normal forcing stress (solid lines; eqn. (12a)) and band-parallel perturbation stress (dashed lines; eqn. (32)). Both are divided by . Calculation details are in Appendix B. (d) and (e) The anisotropy perturbations and divided by , calculated according to eqn. (35).

The growth-rate of porosity perturbations is shown in panel (a) of Figure 6 for and for values of ranging from zero to two. Although does not exclude linearised variations in , comparison with the -curve in Fig. 5b confirms that variations in are ineffectual; with the expected band angle is . For increasing , the growth-rate peak again splits into peaks at low and high angles. It is important to note that the mean value of () is not changed in this exercise — only the amplitude of variations about that mean. Consistent with the numerical results of Figure 1, dynamic variations in the magnitude of anisotropy can sharply reduce band angles, even at moderate for which the static anisotropy model predicts a high band angle ().

Figure 6b breaks the full growth rate into two parts, each associated with one of the terms of equation (15). Dashed lines, representing the shear instability, are almost unaffected by . In contrast, the tensile instability is strongly stabilised with increasing . This stabilisation causes the peak of the full growth rate in panel (a) to split into low- and high-angle peaks. To understand why dynamic anisotropy promotes low band angles, it is therefore sufficient to understand why it stabilises the tensile instability.

The tensile instability is driven by , as discussed in §5.1. This represents the normal stress (tension positive) that arises when viscosity perturbations interact with the base-state strain rate. The detailed form of the forcing stress for the dynamic anisotropy model is given in equation (31). Figure 6c shows that the forcing normal stress varies significantly with (whereas the forcing shear stress, not shown, is almost unaffected by ). The system compliances, which are leading-order quantities, are not affected by dynamic anisotropy. Therefore, it is the variation of that is responsible for stabilisation of the tensile instability under dynamic anisotropy.

To develop a physical understanding of the detailed dependence of on , and (eqn. (31a)), focus attention on , as this is the dominant band angle when . For bands at , solid curves in panel (c) of Figure 6 show that the forcing stress goes from a positive perturbation (in phase with ) to a negative perturbation (anti-phased with ) with increasing — hence the forcing stress in the high-porosity bands goes from tensile to compressive. This change is due to an increase in the magnitude of anisotropy perturbation , shown in panel (d). Since , the deviatoric stress perturbation is entirely due to the band-parallel normal stress perturbation (according to eqn. (32)), which is shown by dashed curves in panel (c). Because , signifies a magnitude reduction of in the high-porosity bands; the largest change occurs for bands at . As sensitivity to deviatoric stress increases, becomes more negative (panel (d)). Negative values of (anti-phased with ) mean high-porosity bands have lower deviatoric stress and weaker anisotropy than the low-porosity, inter-band regions. This is consistent with numerical results in Fig. 1d.

Figure 5d–f shows that increases the normal compliance at angles between zero and 90. A negative perturbation to therefore makes the high-porosity bands in this range of angles less compliant to tensile stress and the low-porosity inter-bands more compliant. Overall, then, the perturbation in anisotropy amplitude tends to cancel the direct effect of the porosity perturbation on the normal compliance, and hence works to stabilise the tensile instability.

The comparison between the static and dynamic anisotropy models developed in this section is summarised in Table 1. These two models are identical at leading order but different at first order. Therefore, stabilisation of the tensile mechanism due to more compressive base-state stress () and destabilisation of the shear mechanism due to shear stress-induced expansion () occur in both the static and the dynamic anisotropy model. These two cases differ, however, in that further stabilisation of the tensile mechanism occurs due to the dynamic variation of anisotropy magnitude (). This effect hardens the band regions and weakens the inter-band regions under dynamic anisotropy. This additional factor ( in Table 1) significantly lowers the band angle.

It is interesting to note that dynamic perturbations to the angle of anisotropy are not an important control on band angle. Figure 6e shows that they are not affected by . More importantly, is always zero for bands orientated at . This indicates that the stabilisation of the tensile instability and the lowering of band angle under dynamic anisotropy cannot be attributed to . In numerical simulations (Fig. 1e), the variations of do not contribute to the lowering of band angle that is observed in Fig 1c, though they are well-explained by the stability analysis at (red dashed line).

5.3 The effect of contiguity increase in the -direction

Until now, we have neglected and focused on the effects of , which quantifies contiguity decrease in the direction of maximum tension. Non-zero represents a weakening in the direction that (i) reduces the magnitude of tensile stress and leads to (ii) shear strain-induced compression and (iii) shear stress-induced extension (Takei and Katz, 2013). We have shown that the tensile mechanism is stabilised around by the first of these and is stabilised around and by the second; we have also shown that the shear mechanism is destabilised by (iii). The parameter quantifies the contiguity increase in the direction of maximum compression. Even if is zero, a non-zero creates viscous anisotropy (see eqn. (2b)), causing the couplings (ii) and (iii). However, Fig. 4 shows that does not cause the weakening (i). It is this weakening, by only, that is responsible for splitting the growth-rate peak in both static and dynamic models. On this basis, we expect the effect of to be small. This is indeed the case: as shown below, alone does not cause a lowering of band angle, but it can affect the lowering by .

Figure 7: The effect of on the growth rate of porosity perturbations for static (left column) and dynamic (right column) anisotropy. (a) Full growth rate from eqn. (21) for and various values of . A single curve for and is also shown. In panel (b), the cases with are decomposed into tensile and shear parts. Line greyscale has the same meaning as in panel (a); there is no curve for . (c) Full growth rate from eqn. (15) for and various values of . A single curve for and is also plotted. In panel (d), cases with are decomposed into tensile and shear parts. (e) The band-parallel, normal stress perturbation .

The left column of Figure 7 illustrates the effect of under static anisotropy. Panel (a) shows that under static anisotropy, is split into high and low angle peaks for any value of when (solid curves), whereas it is peaked at for any value of when ( shown by dash-dot curve). For , increasing causes a modest shift to more compressive at and and a modest increase in the amplitude of shear stress (Fig. 3, right column). Therefore, as Figure 7b shows, causes stabilisation of the tensile instability and destabilisation of shear instability in equal measure. These two effects compensate each other and the solid growth-rate curves in Figure 7a are thus all very similar to that for .

The right column of Figure 7 shows how affects dynamic anisotropy. Panel (c) shows that for in the dynamic anisotropy model, a two-peaked growth rate occurs for and (light gray curve) but does not for and (dash-dot curve). In the former case of non-zero with a double peak, increasing enhances the stabilisation of the tensile mechanism at and deepens the valley between low-angle and high-angle peaks of (Fig. 7c–d). This occurs because is enhanced by the overall strengthening effect of (Fig. 7e). The very low band angles that emerge in the numerical simulation with dynamic anisotropy and (Fig. 1) are therefore a consequence of both the dynamic effect of and the enhancement by .

6 Summary and discussion

We have developed and analysed a model of coupled magma/mantle dynamics with anisotropic viscosity. The anisotropy is controlled by the orientation of principal stresses and the amount of deviatoric stress. The model presented here introduces small modifications on that of Takei and Katz (2013); in particular, the parameter models an increase in contiguity of grains in the direction of maximum compressive stress and the parameter allows for a finer control on the magnitude of anisotropy and its sensitivity to stress (for ). This description of viscous anisotropy is physically consistent with experiments and relatively simple, so its analysis should clarify the mechanics of rocks for which the assumptions hold. Existing experimental data, however, are not enough to quantitatively constrain all parameter values. The parameter studies performed here aim to understand the underlying physics.

It is known from previous theoretical work th