Magmatic channelisation by reactive and shear-driven instabilities at mid-ocean ridges: a combined analysis

It is generally accepted that melt extraction from the mantle at mid-ocean ridges (MORs) is concentrated in narrow regions of elevated melt fraction called channels. Two feedback mechanisms have been proposed to explain why these channels grow by linear instability: shear flow of partially molten mantle and reactive flow of the ascending magma. These two mechanisms have been studied extensively, in isolation from each other, through theory and laboratory experiments as well as field and geophysical observations. Here, we develop a consistent theory that accounts for both proposed mechanisms and allows us to weigh their relative contributions. We show that interaction of the two feedback mechanisms is insignificant and that the total linear growth rate of channels is well-approximated by summing their independent growth rates. Furthermore, we explain how their competition is governed by the orientation of channels with respect to gravity and mantle shear. By itself, analysis of the reaction-infiltration instability predicts the formation of tube-shaped channels. We show that with the addition of even a small amount of extension in the horizontal, the combined instability favours tabular channels, consistent with the observed morphology of dunite bodies in ophiolites. We apply the new theory to MORs by calculating the accumulated growth and rotation of channels along streamlines of the solid flow. We show that reactive flow is the dominant mechanism deep beneath the ridge axis, where the most unstable orientation of high-porosity channels is sub-vertical. Channels are then rotated by the solid flow away from the vertical. The contribution of the shear-driven instability is confined to the margins of the melting region. Within the limitations of our study, the shear-driven feedback is not responsible for significant melt focusing or for shallowly dipping seismic anisotropy [abridged].


INTRODUCTION
At mid-ocean ridges, plate spreading induces upwelling of the mantle, causing decompression melting. This melting occurs inside a volume (the melting region) that extends to a depth and distance of order 100 km from the ridge axis. The rock within the volume is partially molten, consisting of a crystalline solid with liquid melt along the boundaries of solid grains. The melt resides in an interconnected, permeable network of pores, such that it can segregate from the residue and migrate over large distances, driven by buoyancy and pressure gradients. Several lines of evidence suggest that the migration is not spatially uniform, but is rather localised in channels of elevated melt fraction (porosity).
A key line of evidence for channelised transport comes from geological observations of tabular bodies of nearly pure olivine (dunite) in ophiolites, which are otherwise dominantly of olivine+pyroxene lithology (harzburgite). These have been interpreted as the relics of former channels in which focused melt flow has dissolved all pyroxene and replaced it with olivine (Quick 1982;Kelemen 1990;Kelemen et al. 1995aKelemen et al. , 1997Kelemen et al. , 2000Braun & Kelemen 2002). Similar features are observed in laboratory experiments in which Si-undersaturated melts are forced to traverse a porous, olivine+pyroxene matrix (Pec et al. 2015(Pec et al. , 2017. Furthermore, there is well-documented chemical disequilibrium between erupted lavas and the harzburgitic uppermost mantle. This surprising observation can be reconciled by channelized flow; transport through dunite channels isolates rising melts from the surrounding mantle harzburgite (Kelemen et al. 1995a). Moreover, young lavas are typically found to contain isotopes of the uranium-series decay chain with disequilibrium activity ratios . This has been interpreted to indicate rapid melt transport from depth in the mantle, which would be a consequence of the hypothetically channelised flow Elliott & Spiegelman 2014).
Physical models and laboratory experiments point to two possible types of fluid-mechanical instability that can cause localisation of magmatic flow. First, chemical reactions during magmatic ascent can dissolve the host rock due to a solubility gradient, driving a reaction-infiltration instability (Kelemen 1990;Aharonov et al. 1995;Kelemen et al. 1995b). Second, shear of the solid rock coupled with the fact that its viscous strength decreases with melt fraction cause a shear-driven banding instability (Stevenson 1989;Holtzman et al. 2003a). These mechanisms have been extensively studied in isolation, and their outcomes have been invoked separately to explain natural observations.
The key objective of the present paper is to assess the relative contributions of these two instabilities to melt transport within the mid-ocean ridge melting region. In particular, it is to discover the dynamic and parametric conditions under which the shear-driven instability can contribute significantly to the overall pattern of channelised melt transport. It has been widely assumed that the shear-driven instability makes an important contribution to melt flow and mantle dynamics (e.g., Kohlstedt & Holtzman 2009;Kawakatsu et al. 2009;Holtzman & Kendall 2010), but no quantitative assessment that also includes reactive flow has been published. In the context of our simplified analysis, we test and ultimately challenge this assumption. Instead we argue that the reaction-infiltration instability is dominant and that shear-driven instability may be insignificant beneath mid-ocean ridges.
We reviewed research into the reaction-infiltration instability in Rees Jones & Katz (2018) and augment that review briefly here. Recent laboratory experiments at high temperature and pressure showed that highly permeable, cylindrical conduits form due to the reactive flow of Siundersaturated melt (Pec et al. 2015(Pec et al. , 2017(Pec et al. , 2020. Similar cylindrical features also arise from reactive instabilities in a class of porous media called mushy layers (Tait et al. 1992;Worster 1997). The cylindrical geometry of channels in mushy layers is noteworthy because it contrasts with the tabular nature of dunite bodies in ophiolites. Theoretical work has focused on linear stability analysis and twodimensional numerical calculations (Aharonov et al. 1995;Spiegelman et al. 2001;Hewitt 2010;Liang et al. 2010;Hesse et al. 2011;Schiemenz et al. 2011;Szymczak & Ladd 2013Jordan & Hesse 2015). Numerical calculations of reactive flow were extended by Baltzell et al. (2015) to include the effect of mantle shear. However, that study used a constant shear viscosity rather than one that is porosityweakening. Thus it considered the effect of shear in stretching and rotating high-porosity features caused by the reactive instability, but excluded the potential for the shear flow to drive a melt-localising instability.
The shear-driven instability was reviewed by Kohlstedt & Holtzman (2009). The instability was predicted theoretically by Stevenson (1989) before being confirmed experimentally by Holtzman et al. (2003a). It has since been the subject of extensive laboratory experiments (Holtzman & Kohlstedt 2007;King et al. 2011b;Qi et al. 2015) and theoretical study (Spiegelman 2003;Katz et al. 2006;Butler 2009;Alisic et al. 2016). Recent consensus is that it is controlled by viscous anisotropy (Takei & Holtzman 2009a,b,c;Butler 2012;Qi et al. 2015), although the details remain incompletely understood. The emergent patterns, including the associated development of a distinct mode of olivine lattice preferred orientation (Holtzman et al. 2003b) and hence seismic anisotropy (Holtzman & Kendall 2010), have been invoked to explain proxy measurements.
A set of laboratory experiments reported by King et al. (2011c) considered the combined role of reaction and shear in generating melt bands. However, the experimental conditions are so far from the natural system that it is difficult to draw general conclusions from that study.
The natural system of primary interest here is the midocean ridge (MOR). MORs are a fundamental component of plate tectonics, the predominant locus of present terrestrial magmatism, and a context in which melt channelization is inferred from observations, as discussed above. The stability analysis of Butler (2009) suggested that shear-driven porosity bands would form here. However, that study neglected the role of buoyancy and ongoing melting. Vestrum & Butler (2020) considered the effect of both buoyancy and a uniform background melting rate (but excluded reactive melting). They showed that the consequences of uniform melting depend on details of the rheological model. Furthermore, they showed that buoyancy does not affect the magnitude of the shear bands; rather, it causes bands to travel as porosity waves. These studies cannot, however, address the relative importance of reactive flow to melt localisation. The relative importance was tested in terms of its geochemical consequences by Liu & Liang (2019), who considered models with various distributions of both reactive channels and shear-induced bands. However, the distribution and character of localised porosity features in Liu & Liang (2019) were prescribed, rather than emerging dynamically.
The broad goal of this paper is to develop a theoretical understanding of the combined dynamics of reactive and shear-driven instabilities in the partially molten upper mantle. This allows us to assess their relevance for magmatic flow localisation at MORs, at length-and time-scales that are necessarily very different from those of the laboratory. Insights gained here may have wider implications for magmatism.
In section 2, we develop a theoretical method to determine the linear growth rate of perturbations in the form of alternating bands of higher and lower porosity. The theory allows for the perturbations to evolve by both reactive and shear mechanisms simultaneously. We show that this theory can reproduce previous estimates of their growth rate in both the reaction-only and shear-only limits.
In section 3, we apply this theory to the idealized scenario of an infinite, partially molten material with a uniform background magmatic flow, a linear solubility gradient aligned with gravity, and a linear shear of the solid matrix. We identify a parameter that describes the relative importance of reaction and shear. In this idealised context we show that there are two distinguished directions that control the orientation of high-porosity pathways: one in the direction of gravity (the vertical) over which there is a chemical solubility gradient, and the other in the direction of the maximum tensile deviatoric stress. We calculate the growth of the instability and show that the optimal orientation for growth of porosity perturbations depends on these directions and also the ratio of shear-to-reactive growth rates. Indeed, within a plane normal to gravity, it is shear that breaks the horizontal isotropy. We show that the most unstable, fastest growing features are tabular bodies extending horizontally in the direction in which there is no component of the shear flow. This is important because dunites are observed to be tabular features, rather than the cylindrical conduits favoured by the pure reaction-infiltration instability. In this section we also discuss the role of compaction, chemical advection and diffusion, all of which play a role in determining the wavelength-dependence of the growth rate.
In section 4 we present a methodology and in section 5 we present results from the first combined assessment of both reactive and shear-driven instability at mid-ocean ridges. We build on the approach of Gebhardt & Butler (2016), allowing the amplitude and orientation of perturbations to evolve along streamlines of the solid flow, based on their local growth rate. But whereas Gebhardt & Butler (2016) considered only shear-driven growth, in our case, the growth rate is a consequence of both reaction and shear. We show that their contributions are not spatially uniform and predict how they vary along corner-flow streamlines. Generally, we find that the initial growth of porosity bands is dominated by reaction. Within the space of the melting region and the plausible ranges of control parameters, our results indicate that melt channels are sub-vertical features that undergo some rotation by the shear flow. Shear-driven instability is predicted to contribute at the margins of the melting region.
In section 6, we review the implications of these findings for understanding the relative importance of reaction and shear in driving channelized melt extraction from the mantle. We discuss implications for interpreting the origins of tabular dunites embedded in a harzburgitic upper mantle. We also discuss melt focusing towards the ridge axis, as well as seismic anisotropy that has been attributed to aligned, melt-rich bands.

Equations governing two-phase flow
The partially molten upper mantle can be modelled as a region of two-phase flow. The continuum model that we use, developed by McKenzie (1984), is based on averaging all the quantities of interest across a control volume containing both solid and liquid phase. Table 1 lists the state variables  and table 2 lists the physical properties that we define in this section. Our formulation makes a Boussinesq approximation: the solid density ρs and liquid density ρ l are taken to be constants, and their difference ∆ρ = ρs − ρ l is neglected everywhere, except in so far as it drives segregation of the liquid by buoyancy. Mass conservation in the liquid phase is given by where t is time, φ is the porosity, v l is the liquid velocity and Γ is the volumetric melting rate. Similarly, mass conservation in the solid phase is given by where vs is the solid velocity. The compaction rate is defined as C = ∇ · vs. It is convenient to sum equations (1) and (2) to form an equation for the compaction rate: where the Darcy (melt segregation) velocity is defined as Then equation (2) can be rewritten In this formulation, porosity, moving with the solid phase, evolves due to melting rate Γ and compaction C. Note that C > 0 represents a decompaction of the mantle matrix, i.e., an increase in the porosity. Equations (3) and (5) define two-phase mass conservation in our system. Next, two-phase momentum conservation can be written using a 'Stokes-Darcy' formulation. The Darcy velocity is driven by gradients in the liquid pressure P l and buoyancy (g is gravity). K is the liquid mobility, which is defined as the permeability divided by the liquid viscosity. We will often refer to K as the permeability, since we assume the liquid viscosity is a constant so variation in K comes from variation in permeability. The Stokes part of the system can be written where η is the shear viscosity (Newtonian), ζ is the bulk viscosity, ρ = φρ l + (1 − φ)ρs = ρ l + (1 − φ)∆ρ is the bulk density and is the deviatoric strain rate tensor (I is the identity tensor, T is the transpose operator). We discuss constitutive laws and the dependence of material properties on porosity in the next section. For now, we note that the shear-driven mode of instability relies on the fact that η decreases with φ. This motivates expanding out the derivatives involving η. We also apply a Helmholtz decomposition of the solid velocity into a shearing (incompressible) part and a compacting (compressible, but irrotational) part where U is a scalar potential that can be related to C by the relationship C = ∇ 2 U . Then equation (7) becomes: The terms involving pressure, η∇ 2 u and gravity can be recognised as the usual, single-phase form of Stokes law with a constant Newtonian viscosity. We can eliminate various terms by taking the curl of equation (10) to form a type of vorticity equation We can also substitute equation (10) into equation (6) to obtain vD = −K 2Ds · ∇η + η∇ 2 u + 4 3 η∇C + ∇(ζC) + (1 − φ)∆ρg .
(12) The melting rate Γ is determined by chemical disequilibrium. We describe the model that we and others have used in more detail in Rees Jones & Katz (2018); the original model is due to Aharonov et al. (1995). The crucial ingredients are that the melting rate is linearly proportional to the chemical disequilibrium where R is the reaction rate constant and X is the chemical undersaturation. That R is constant is a reasonable assumption during the initial growth of the instability; later R will decrease as the concentration of the soluble component in the solid phase decreases. Equation (13) is a first order, linear, kinetic reaction rate equation. Then X is governed by an advection-diffusion-reaction equation, which can be written in the form where β is the constant gradient of the equilibrium chemical concentration, sometimes called the solubility gradient, which we assume is orientated in the vertical direction. Under these assumptions, the liquid concentration can be written in terms of the undersaturation as c l = βz − X, appearing in the second term in equation (14). Also, α represents the supersaturation of the reactively produced melts that drive the system back towards equilibrium. Within equation (14), α is an inverse reactivity, because the reactive melt rate Γ is inversely proportional to α. Finally, DX is the diffusivity of chemical species in the liquid phase, scaled by the porosity φ. Diffusion in the liquid is much faster than that in the solid, so diffusivity of chemical species in the solid phase is neglected. All the approximations made here are described in more detail in Rees Jones & Katz (2018). For now we focus on the equilibrium concentration gradient because this drives the reactive instability. In a more complete description, the equilibrium chemical concentration depends on pressure (Kelemen et al. 1995b;Longhi 2002). Thus this simple formulation with constant β combines an assumption that the pressure is dominantly lithostatic, i.e., proportional to vertical position z, and an assumption that the dependence of equilibrium chemical concentration on pressure is linear. The latter assumption could be relaxed straightforwardly by allowing β to vary with z. In appendix A, we consider the full pressuredependence of the equilibrium concentration gradient and hence the reactive melting rate. The shear-driven instability creates a pressure gradient that can, in turn, feed back on the reactive instability through the pressure-dependent reactive melt rate. So appendix A also considers a mode of coupling between the two types of instability. We find the effects of pressure-dependence on the linear growth rate of the reactive instability are relatively small, because ∆ρ/ρs 1, which means that the pressure is dominantly lithostatic. However, it is possible that the pressure-dependence may have a larger effect on the nonlinear development of channels. Spiegelman et al. (2001) suggested that lateral pressure gradients that focus flow towards a mid-ocean-ridge (MOR) axis could promote the development of a diagonal solubility gradient that would further enhance convergence of flow towards the ridge axis. Nonlinear calculations extending those of Spiegelman et al. (2001) could assess the significance of this proposed mechanism at MORs.

Equations governing linear growth of instabilities
Our first goal is to determine the local growth rate of instabilities about a uniform background state with a linear incompressible shear flow u and a vertical Darcy flow wD0z, wherez is a unit vector in the z-direction. We write all the fields as an expansion into a base state (subscript 0) and a perturbation ( ) Table 3 summarizes the notation introduced in this section. A crucial ingredient of the reactive-and shear-driven instabilities is that the constitutive laws (material properties) depend on porosity. We can linearize all the constitutive laws by expanding them in a Taylor series in φ truncated at first order where all the primed variables are proportional to φ . For example A common choice of constitutive law for permeability, and hence melt mobility assuming melt viscosity is constant, is K = K * φ n , where K * is a material property and n is an exponent (e.g. von Bargen & Waff 1986). Then Likewise for shear viscosity, a commonly used law is η = η0 exp [−λ * (φ − φ0)], where λ * is a material property (e.g. Kelemen et al. 1997;Mei et al. 2002). Then The same methodology can be applied to any form of constitutive law that has a continuous first derivative. We do not state any particular constitutive law for the bulk viscosity at this stage, because ζ does not affect the linear analysis that follows. We discuss this issue further in section 4.1.2.

Background state
We are interested in the development of perturbations to an initially uniform porosity field with a uniform, background upward flow of magma and a linear, incompressible mantle shear flow. This state is called the background or base state. A somewhat subtle issue is that a uniform background state is only an approximate solution of equations (16)(17)(18)(19). In Rees Jones & Katz (2018), we showed that this approximation holds provided we are dealing with length scales much smaller than αβ −1 ≈ 500 km. Physically, the system is only weakly reactive. In this approximation, the background compaction rate is negligible, so C0 = 0. This approximation is equivalent to saying that the background melting rate (which balances the background compaction rate) is negligible, an issue we discuss later in the context of mid-ocean ridges (section 4.1.2).
Equation (18) gives an expression for the background Darcy velocity that can also be written in terms of the liquid velocity w0 using wD0 = φ0w0. In particular, Any linear incompressible shear flow satisfies equation (17), because ∇η = 0 and ∇C = 0. We define which is the symmetric velocity gradient tensor associated with the background shear flow (note this is a constant and trace-free tensor, because the background shear flow is linear and incompressible). A mid-ocean ridge far from any offsets is roughly two-dimensional, so if we choose the x-axis to be in the direction of plate spreading, the shear is in the x-z plane. Then a general expression for the symmetric velocity gradient tensor of this type of flow is whereγ0 is the background strain rate and θe is the angle of maximum extension (θe = 0 corresponds to extension in the x-direction and θe = π/2 corresponds to extension in the z-direction). By this choice of co-ordinates, there is no extension or contraction in the y-direction. Equivalently, ψe = 0 where ψ is the azimuthal angle (figure 1).

Linear equations
We now substitute the linear decompositions from equations (20-27) into the governing equations (16-19) and neglect any quadratic terms. The vorticity equation (17) becomes Here, we exploited the fact that the background compaction rate is negligible, which means that the effect of variations in < l a t e x i t s h a 1 _ b a s e 6 4 = " K r n + Q U 7 C 6 b g f O l Y y Z E L 4 L q d k g j k = " > A A A C Y H i c b V F N a x s x E J W 3 X 4 7 7 k b i 9 t Z e l p t C T 2 T U J 7 T G 0 U H p 0 o X Y C 3 s V o t b O x i K Q V 0 m w T I / Q X e m 3 / W q / 9 J d X a S 3 D s D g j e v D c P j Z 4 K L b j F J P n T i x 4 8 f P T 4 S f 9 o 8 P T Z 8 x f H J 8 O X c 1 s 3 h s G M 1 a I 2 l w W 1 I L i C G X I U c K k N U F k I u C i u P 7 f 6 x Q 8 w l t f q O 6 4 1 5 J J e K V 5 x R r G l M m 3 5 8 m S U j J N N x Y c g 7 c C I d D V d D n t f s r J m j Q S F T F B r F 2 m i M X f U I G c C / C B r L G j K r u k V L A J U V I L N 3 W Z Z H 7 8 L T B l X t Q l H Y b x h d x 2 O S m v X s g i T k u L K 7 m s t + T 9 t 0 W D 1 M X d c 6 Q Z B s e 1 F V S N i r O P 2 5 X H J D T A U 6 w A o M z z s G r M V N Z R h y G c w y B T c s F p K q k q X l b y q / G K S u 6 w K A y 5 D u E V X e j d K v d / p J t 7 7 + 0 a 9 7 9 R t K l R 0 z r t u 4 w z J p / s 5 H 4 L 5 Z J y e j s + + n Y 7 O P 3 V / 0 C d v y F v y n q T k A z k n X 8 m U z A g j K / K T / C K / e 3 + j f n Q c D b e j U a / z v C L 3 K n r 9 D w P J u d k = < / l a t e x i t > bulk viscosity turn out to be consequently small (we discuss the effect of variation in bulk viscosity in Rees Jones & Katz (2018) so do not elaborate further here).
Equation (28) motivates the introduction of a scalar potential ψ defined by which (since ∇ · u = 0) must satisfy a Poisson equation Then we introduce a scaled, dimensionless potential where the implication follows using equations (24) and (27). The remaining equations can be rewritten by substituting the constitutive laws and eliminating the Darcy velocity. An expression for the perturbed Darcy velocity is given in appendix A in equation (A.14). The remaining equations become where a subscript z is used to denote a partial derivative with respect to z. Two important parameters emerge: where δ is the compaction length (McKenzie 1984) and Λ controls how much shear-driven compaction arises from a porosity perturbation. The ratio Λ/δ 2 determines the growth rate of the shear-driven instability (section 2.3). This ratio was first identified by Stevenson (1989), in which it is written σm.

Normal modes
We now look for normal mode solutions of the form where σ is the growth rate of the instability, k = [kx, ky, kz] is the wavevector, x = [x, y, z] is the position vector, and the prefactors are constants. We define k as the magnitude of k, so k 2 = k 2 x + k 2 y + k 2 z . The co-ordinate system is shown in figure 1.
We define the function such that equation (31) becomesψ =φG. The use of a factor k −2 in the definition ensures that G is independent of the wavenumber (magnitude of the wavevector) and only depends on the direction. G has a maximum value of +1 when ky = 0 and k is in the direction of maximum extension. G has a minimum value of −1 when ky = 0 and k is perpendicular to the direction of maximum extension. Then equations (32-34) become where is an extended version of the inverse reactivity of the system, which is augmented by advection and diffusion of the undersaturated chemical species.
Finally, we substitute equation (40) and then equation (39) into equation (38) to obtain an expression for the combined growth rate of shear and reactive instabilities In subsequent sections, we explore the nature of this equation in detail and discuss the physical significance of the terms that appear in it. First, we relate it to previous studies of the reactive-and shear-driven instabilities in isolation.

Shear-driven instabilities
The reactive part of the growth rate can be eliminated by setting β = 0. This ensures thatΓ = 0, so the reactive part of the contribution to porosity change is eliminated and we are left with the part coming from shear. Then equation (42) becomes The term involving w0 arising from buoyancy-driven melt flow is typically neglected because buoyancy is unimportant in laboratory experiments that impose a rapid shear (Spiegelman 2003). In any case, it only affects the imaginary part of the growth rate, giving rise to compaction waves. The dependence on the compaction length δ in equation (43) means that when δk 1 (the wavelength of the instability is much greater than the compaction length), then the growth rate approaches zero. Conversely, when δk 1, the real part of the growth rate can be approximated This motivates us to define which is the dimensional growth rate associated with shear (Stevenson 1989). Note that the real part of the growth rate has a distinguished direction, namely the direction of extension in the x-z plane. The imaginary part comes from vertical background magma flow, which drives the perturbations upward as waves. So the vertical is an additional distinguished direction in this case. The growth rate and angular dependence through equation (37) are known from previous studies (e.g., Spiegelman 2003; Katz et al. 2006).

Reaction-driven instabilities
The shear part of the flow can be eliminated by taking Λ = 0. Then equation (42) becomes This can be shown to be equivalent to Rees Jones & Katz (2018) under the same assumption (δk 1) mentioned previously. For now, note that if we additionally make the assumption α ≈ α (valid when the reaction rate R is very fast), then The growth rate has a cylindrical symmetry about the vertical direction. The maximum growth rate occurs when kz = 0, since k 2 x + k 2 y = k 2 − k 2 z . Thus the channels formed by the reaction-infiltration instability are vertical (Rees Jones & Katz 2018). Moreover, we show in appendix A that the preferred vertical orientation holds even when we consider the full pressure-dependence of the solubility gradient (at least in the context of linearised analysis).
In this section we have shown that the full dispersion relation (42) includes, as special cases, the results of previous studies on both the shear-driven instability and the reactioninfiltration instability.

RESULTS: LOCAL ANALYSIS OF COMBINED INSTABILITY
In this section we analyse controls on the growth rate σ that are relevant to the case of an infinite domain with a uniform base state.

Equilibrium dynamics at large compaction length
The simplest version of the instability involving both reaction and shear can be illustrated by considering an important limit of equation (42). In particular, if the reaction rate is very fast, the system is driven to equilibrium and α ≈ α. If also we consider the short-wavelength or large-compactionlength limit discussed earlier (δk 1), then It is important to note that all the terms involving the wavevector are independent of its magnitude; they depend only on its direction. This is a feature of the particular limit considered that (by design) neglects the role of advection, diffusion and compaction in affecting the wavelength. We consider these controls later.
In this particular limit, the behaviour is controlled by the ratio of the growth rate of shear-driven to reactiondriven instabilities, which we can write where we grouped together material parameters separately to the ratioγ0/w0. The former group might be expected to be roughly constant, provided the bulk-to-shear viscosity ratio is constant, whereas the latter will vary spatially at a mid-ocean ridge.

Three-dimensional effects and the orientation of porosity bands
The reactive mode of instability has a cylindrical symmetry (there is no difference between the x and y direction; the only distinguished direction is the vertical z). Numerical calculations show that the instability leads to the formation of cylindrical, high-porosity conduits (M. Spiegelman, unpublished work), in accordance with laboratory experiments (Pec et al. 2015(Pec et al. , 2017. However, the shear-driven mode of instability leads to the formation of high-porosity sheets, the orientation of which depends on the direction θe of maximum rate of extension. The coordinate system is aligned such that shear is in the x-z plane, with the direction of extension having an azimuthal angle ψe = 0 to that plane. Hence porosity sheets extend parallel to the y-direction. In cross-section on the x-z plane, the sheets appear as highporosity bands (which is how laboratory experiments are typically presented). Here, we investigate the combined effect of the reactive and shear mechanisms in light of the fact that that this combination has two distinguished directions: the orientation of the shear flow and the direction of the solubility gradient (vertical). We work in terms of modified spherical polar coordinates (θ, ψ), shown in figure 1a, where −π/2 ≤ θ ≤ π/2 is the inclination angle and 0 ≤ ψ < 2π is the azimuthal angle. This definition of θ is consistent with that of θe (figure 1b). However, it is not the same as the angle of the porosity bands as it was defined in, for example, Spiegelman (2003). Rather, as shown in figure 1a, the porosity bands are normal to the wavevector that makes an angle θ to the x-y plane. In the particular case that we restrict attention to the x-z plane (ψ = 0), the angle of the porosity bands according to the definition of Spiegelman (2003) is π/2 − θ.
(53) For the particular case that ψ = 0, since the growth rate only depends on the wavevector orientation relative to the direction of extension. Thus G = 1 when the wavevector is parallel to the direction of extension and G = −1 when they are perpendicular. Figure 2 shows the growth rate real(σ) from equation (52) as a function of θ (x-axis) and ψ (y-axis). This function has four stationary points; we consider each in turn. (i) There is one stationary point at θ = ψ = π/2, which is always a saddle point. This is a mode with a wavevector purely in the z-direction. Therefore it represents a tabular perturbation aligned with the x-y plane. (ii) There is another stationary point when θ = 0 and ψ = π/2. This is a mode with a wavevector purely in the y-direction, that is, geometrically, a tabular feature in the x-z plane. The nature of this stationary point depends on the magnitude of the shear. When S exceeds a critical value S > Sc ≡ − cos 2θe, it is a saddle. However, when S < Sc, this stationary point is a maximum. (iii/iv ) There are two more stationary points with ψ = 0 (so ky = 0) and tan 2θ = 2S sin 2θe 1 + 2S cos 2θe .
(55) Figure 3 shows these modes, which are tabular features that have a wavevector in the x-z plane at an angle θ; they extend in the y-direction. Equation (55) has two roots, one in the domain −π/2 ≤ θ < 0 and another in the domain 0 ≤ θ < π/2. The former root is associated with contractional shear stress and is always a minimum. The latter root is associated with extensional shear stress and is a maximum when S exceeds the aforementioned critical value S > Sc ≡ − cos 2θe and a saddle when S < Sc.
We expect the most unstable mode (i.e., the one that grows fastest) to be dominantly expressed in a full solution of the governing equations and hence we analyse what selects this mode. Figure 4 shows that there is a transition for a horizontally isotropic (no difference between x and y directions) to an anisotropic growth rate as the shear rate S increases, moving from left to right in the figure. However, the transition depends on the angle of maximum extension in the shear flow (from top to bottom in the figure). If cos 2θe ≥ 0, which corresponds to the rows with θe = 0, π/8, π/4, the transition is immediately to a state where ky = 0. Geometrically, when cos 2θe ≥ 0, the angle of maximum extension is within π/4 of the horizontal, which means that the horizontal direction is extensional and the vertical direction is contractional, as shown by the mantle flow plots in the left-most column of figure 4. This leads to tabular features that are orientated vertically for small S, because reaction dominates the instability and promotes vertical features. The tabular features approach an orientation perpendicular to the angle of maximum extension as S increases, because this orientation is favoured by the sheardriven instability. However, if cos 2θe < 0, which corresponds to the rows with θe = 3π/8, π/2, the situation is reversed: the x-direction is contractional and the vertical direction is extensional, as shown by the mantle flow plots in the leftmost column of figure 4.. Then for small S, the most unstable wavevector orientation has kx = 0, kz = 0 (i.e., a tabular feature aligned with the x-z plane) due to the combined effect of reaction and shear. As in previous cases, it is reaction that promotes vertical features (kz = 0). However, in contrast to previous cases, shear-driven contraction suppresses wavevectors in the x-direction so kx = 0. As S increases through the critical value defined above Sc ≡ − cos 2θe where the shear-driven instability dominates over the reaction-driven instability, the state switches to ky = 0 with orientation approaching the angle of maximum extension for large shear, as before.
Mid-ocean ridges (MORs), which we will turn to in the second half of the paper, are a horizontally extensional environment, corresponding to the upper rows of figure 4. So the crucial implication of this figure for MORs is that even an extremely small shear S > 0 is sufficient to break the horizontal isotropy (symmetry under any coordinate rota-  (52) in the range −π/2 < θ < π/2 (horizontal direction) and 0 < ψ < π (vertical direction). These axis limits are shown in the labelled diagram on the bottom row. The colour scale shows the normalized growth rate relative to the maximum (growth is red, no growth is white, decay is blue). Crosses represent saddle points, squares represent maxima, circles represent minima. Icons on the left-most column indicate the orientation of the background solid shear flow. tion about the z axis) that occurs when only the reactive instability operates. So rather than the expecting the tubeshaped channels that arise from the pure reactive instability (S = 0), the most linearly unstable feature with shear S > 0 are tabular bands. This is important, because dunite channels (interpreted as relics of high porosity channels) have a tabular morphology, as discussed in section 1.

Effect of the compaction length
We now relax the assumption that perturbation wavelength is much shorter than the compaction length. Then equa-tion (52) generalizes to This has the same angular dependence as in the smallwavelength limit kδ 1, so all the conclusions of section 3.2 still hold, including the optimal wavevector orientation. Figure 5 shows how the results depend on the compaction length. It plots the normalised growth rate as a function of kδ for perturbations orientated with ψ = 0 and θ chosen to maximize the growth rate. As in figure 3, there are two regimes depending on whether the z-direction is in extension or contraction. The upper panel shows the case of θe = π/8, which has contraction in the z-direction. Here, Figure 3. Cross section of the normalized growth rate from figure 2 for modes restricted to the x-z plane (i.e., ky = 0) for θe = π/8 (upper) and θe = 3π/8 (lower). Note that the normalization is relative to the maximum growth rate for this restricted set of modes, a distinction only relevant where the stationary point is a saddle point (crosses) rather than a maximum (squares).
the contributions to the growth rate from reaction and shear vary oppositely with the compaction length. The reactive contribution, which is the first term on the right-hand-side of equation (56), decreases slightly with kδ. The amount of decrease depends indirectly on S, because for larger values of S, the angle θ that maximises the overall growth rate increases from 0 towards θe. This decreases cos 2 θ in the first term of (56). By contrast, the shear-driven contribution, which is the second term on the right-hand-side of equation (56), increases significantly with kδ, starting from zero when kδ 1. Thus the total growth rate is dominated by reaction when kδ is small, but there can be a cross-over as kδ increases. This cross-over depends on S. Indeed when S is very small (e.g., the solid curve for S = 0.2) reaction always dominates. Physically, the shear-driven instability is suppressed when the wavelength is comparable to the compaction length (section 2.3), while the reaction-driven instability is not directly affected by compaction in the equilibrium limit ( α = α). Disequilibrium effects introduce a direct dependence on the compaction length (Rees Jones & Katz 2018).
The lower panel shows the case of θe = 3π/8, which has contraction in the x-direction. The behaviour is different when S < Sc = − cos 2θe, as illustrated by the solid curves for S = 0.2. Unlike the upper panel, the shear-driven contribution decreases with kδ, starting from zero when kδ 1 and decreasing towards a negative value when kδ 1. This is because of the contraction in the direction of the optimal wavevector, consistent with figure 3. Thus the total growth rate decreases with compaction length in this case. However, for S ≥ Sc, there is extension in the direction of the optimal wavevector, and the system behaves in the same way as in the upper panel.

Disequilibrium effects
So far, we made the assumption that the reaction rate was extremely fast such that α ≈ α. We now consider the role of disequilibrium effects caused by advective or diffusive chemical transport. We re-write equation (41) The two dimensionless parameter groups are Damköhler numbers that express the reaction rate relative to advective and diffusive transport, respectively. These are typically extremely large and hence the system is typically close to equilibrium (Aharonov et al. 1995;Spiegelman et al. 2001;Rees Jones & Katz 2018). Physically, advection of undersaturation is negligible when Daw → ∞ and diffusion is negligible when DaD → ∞.
We first restrict attention to the role of advection by taking the limit DaD → ∞ so We substitute this expression into equation (42) and find where We make the approximation 1 + Da −2 w sin 2 θk 2 δ 2 ≈ 1 and also take the limit δk 1. Then equation (60) simplifies to This expression can be compared to equation (52) and shows that the advection of undersaturated melts in the presence of shear causes an additional contribution to the growth rate. The contribution is proportional to the shear growth rate, but much smaller, because βφ0w0/α 2 R 1, so Sr S. Physically, this small contribution comes from the shearinduced perturbation melt flow advecting against the background equilibrium concentration gradient combined with the background melt flow advecting the perturbed undersaturation [see left-hand-side of equation (19)]. Although this contribution is small, it is physically interesting because it arises from the coupling of the shear and reactive instability, rather than from the sum of their separate rates.
We next turn attention to the role of diffusion by taking the opposite limit Daw → ∞ so  (52) projected on to the kx-ky plane, as indicated on the circular icon on the bottom. The growth rate is calculated as a function of (θ, ψ) and then plotted against [kx(θ, ψ), ky(θ, ψ)]. The result in independent of the magnitude of the wavevector. The colour scale shows the normalized growth rate relative to the maximum (red is growth; blue is decay; white is neutral). Squares represent maxima. Squares are not shown for the column S = 0; in this case, the maximum growth rate is achieved on the circle k 2 x + k 2 y = k 2 . Icons on the left-most column indicate the orientation of the background solid shear flow.
In this case there is no new contribution to the growth rate; the only change is that the reactive contribution to growth is slightly reduced, especially at high wavenumber (because diffusion acts at small length scales), as discussed in Rees Jones & Katz (2018).
In the next section we embed the above considerations of perturbation growth into the background, large-scale flow beneath a mid-ocean ridge. We neglect both disequilibrium effects and also the role of compaction length (equivalently, the role of the wavenumber). Figure 5 indicates that the growth rate depends weakly on wavenumber provided (kδ) 2 100, or equivalently provided the wavelength is less than about δ × 10/2π ≈ 1.6 km, for a compaction length δ ≈ 1 km. Given that observed dunite channels are smaller than this (Braun & Kelemen 2002), we make the simplifying assumption that kδ 1. Figure 5. Dependence on normalized compaction length kδ. As in figure 3, the upper panel shows θe = π/8 and the lower panel shows θe = 3π/8. In both panels, we plot equation (56) for ψ = 0 and θ chosen to maximize the growth rate according to equation (55).

METHODS: GROWTH OF INSTABILITIES AT MORs
At mid-ocean ridges (MORs), plate spreading drives a circulation of the upper mantle. This is a viscous shear flow and so could promote the formation of shear-driven porosity bands. The upwelling of the mantle beneath the ridge causes decompression melting in a roughly triangular region down to a depth of around 80 km (Langmuir et al. 1992). The resulting melt upwells due to its buoyancy. Over this depth (and hence pressure) range, there is a gradient in the equilibrium chemistry of magma that can drive reactive melting (Kelemen et al. 1992(Kelemen et al. , 1995aLonghi 2002). Thus mid-ocean ridges are a geological setting where both shear-driven and reaction-driven porosity localisation have the potential to occur. In this section, we estimate the relative and combined contributions of these mechanisms. We discuss the predicted orientation of the localized features that result from perturbation growth.

Background state of the mid-ocean ridge
Our approach combines models of several aspects of midocean ridges with the calculations of the growth rate of porosity bands made above. First we estimate the background state, that is, the behaviour of the system in the absence of any small-scale localization of porosity. Our esti-mates focus on the partially molten region beneath the lithosphere, shown in figure 6. For present purposes, the lithosphere is defined as a rigid plate, moving uniformly away from the ridge axis at a speed U0. In section 4.1.1, we calculate the passive flow of the partially molten mantle as the response to motion of the lithosphere. In section 4.1.2, we calculate the background magma flow and porosity by dividing the partially-molten region into a series of melting columns in which magma rises vertically. We assume that each of these melting columns terminates at the base of the lithosphere. This simple formulation for the background-state magma flow precludes capture of lateral flows associated with gradients in compaction pressure (e.g., melting-ratepressure focusing (Turner et al. 2017;Sim et al. 2020), focusing within a decompaction channel immediately beneath the lithosphere (Sparks & Parmentier 1991;Spiegelman 1993;Ghods & Arkani-Hamed 2000)). However, if melt transport throughout most of the melting region (aside from a narrow region near the base of the lithosphere) is buoyancy-driven and hence nearly vertical, our representation is valid, if approximate. We return to this issue in section 6.3.

Mantle flow and strain rate
We estimate the background flow by treating the mantle as a uniform, isoviscous material that flows in response to prescribed plate motion. This is sometimes called a passive or kinematic model. In treating the mantle as isoviscous, we are neglecting the variation of shear viscosity with the background porosity. Based on the parameters estimated subsequently (λ * = 26, φ0 2 × 10 −3 ), this is a reasonable approximation since e −26×2×10 −3 ≈ 0.95, i.e., only about a 5% reduction in viscosity. Equations governing the flow of an isoviscous material in a triangular region are given by Batchelor (1967) and applied to a mid-ocean ridge melting region by Spiegelman & McKenzie (1987). Here, we summarize these results. Figure 6(a) shows the co-ordinate system, centred on the ridge axis, with x the horizontal distance from the ridge and z the vertical distance from the surface, measured upwards, such that z = −H is the bottom of the melting region. The corner-flow solution is most naturally expressed in polar co-ordinates (r, ϑ), where r is the distance to the origin and ϑ is the angle to the downward vertical (so ϑ = 0 is straight down). The partially molten region is triangular and extends to ϑ = ±ϑ1, where π/2 − ϑ1 is the dip of the bottom of the lithosphere.
There is a separable solution in this geometry where the radial flow ur and tangential flow u ϑ can be written The function Θ(ϑ) that satisfies the relevant boundary conditions (that the flow is symmetric about the ridge axis and uniformly translating at speed U0 in the lithosphere) is where C = ϑ1 − sin ϑ1 cos ϑ1 is a constant that depends only on the geometry of the melting region. We can convert from polar to scaled Cartesian coordinatesz = −r cos ϑ,x =r sin ϑ, where all distances are scaled by H, e.g.,r = r/H. Then the scaled solid velocitỹ u = u/U0 has components The scaled velocity gradient tensor is which has a symmetric part Then, by comparison with equation (27), we find that the dimensional strain rate iṡ and the direction of maximum extension is Thus when x < −z, the angle of maximum extension is below the horizontal (θe < 0) and vice-versa when x > −z.
The strain rate is zero beneath the ridge axis (x = 0) and increases in approach to the base of the lithosphere (as x increases or r decreases). The term involving sign(ϑ) in equation (71) arises becauseγ > 0 by definition, leading to the the modulus operator in equation (70). In cylindrical polar co-ordinates, the only non-zero components of the symmetric rate-of-strain tensor are D rϑ = D ϑr =γ (which can be seen from the separable solution in terms of Θ(ϑ)). Figure 7(a) shows this background strain rate. The other information that we need for the magma flow calculations in the following section is the mantle upwelling speed W b at the base of the melting column (wherez = −1 andr 2 = 1 +x 2 ). Here, (1 +x 2 ) −1 − cos 2 ϑ1 C .
As expected, this equation says that the upwelling speed decreases with distance from the ridge axis.

Magma flow and porosity
We estimate the background magma flow and porosity by dividing the melting region into a series of one-dimensional melting columns. The methodology for calculating the flow in a melting column was originally devised by Ribe (1985). The presentation and notation is based on Rees Jones & Rudge (2020), which develops a revised estimate of magma velocities based on the magmatic response to the deglaciation of Iceland (Maclennan et al. 2002). We start from the equations of two-phase flow presented in section 2.1. Making the same simplifications and taking the special case of a one-dimensional flow at steady state, equations (15, 16) imply that where w is the vertical component of the magma velocity. Then we take Γ = Γ0, where Γ0 is the rate of decompression melting, which we assume to be uniform within the melting column. This is a reasonable approximation in the present context, but it neglects the drop in productivity when a mineral phase is exhausted from the residue (Hirschmann et al. 1999). Moreover, it does not apply for volatile-driven melting at depths beneath the dry (i.e., volatile-free) solidus, or at the transition to dry melting (Keller & Katz 2016).
As an aside, note that in our earlier development of the method in section 2 we could have used an alternative melting-rate parameterization. In place of equation (19), we could have used in which case Γ0 would have entered the base state calculations but would not have affected the linear equations governing the perturbations. The most important difference is that the background compaction rate would no longer have been negligible, which would have complicated the analysis. Physically, the background compaction rate is a stabilizing influence on the reaction-infiltration instability, the magnitude of which is very sensitive to the dependence of bulk viscosity on porosity (  The parameterization of equation (74) also neglects the disequilibrium effects (section 3.4). Thus we use equation (19) to calculate the melting rate throughout this study, which maintains consistency with previous studies of the reactioninfiltration instability (Aharonov et al. 1995;Rees Jones & Katz 2018). Equation (73) can be integrated, with the constant of integration chosen such that φw = 0 at the base of the melting column z = −H, in which we made use of the fact that Γ0H = FmaxW b , where Fmax is the maximum degree of melting and W b is the mantle upwelling velocity at the base of the column. Then we take equation (18) and assume that background magmatic segregation is entirely driven by buoyancy to obtain φw = Q0φ n , Q0 = K * ∆ρg, where K * is the prefactor in the permeability-porosity relationship K = K * φ n , and ∆ρg is the buoyancy associated with the density difference between solid and liquid phases (see section 2.2). We combine equations (75) and (76) to We combine this last expression with equation (72) to obtain an expression for the background rate of magmatic upwelling relative to the plate half-spreading rate, (79) Here, we relabelled w as w0, using a subscript to refer to a background property, consistent with the convention introduced in section 2.2. Figure 7(b) plots the ratio w0/U0, and figure 7(c) shows the corresponding background porosity, denoted φ0. These plots are sensitive to the efficiency of melt extraction Q0. If melt extraction is more efficient (Q0 is higher), then the relative melt velocity w0/U0 is increased, reflecting the proportionality to Q 1/n 0 in equation (79). The steady-state porosity φ0 is a balance between melt production and melt extraction. Thus if melt extraction is more efficient (Q0 is higher), then φ0 is decreased, reflecting the proportionality to Q −1/n 0 in equation (77).

Choice of parameters
Our calculations use a reference set of parameters to illustrate the typical behaviour of the model. This includes a prescribed dip of the base of the lithosphere of π/5 (36 • ), as an intermediate case between a very steep boundary (π/4, or 45 • ) and a much shallower case. Later, in section 5.4, we consider the behaviour when π/12 (15 • ), which might be appropriate if the lithosphere is interpreted to be the cold thermal boundary, which thickens gradually at intermediate to fast spreading rates.
Reference material parameters are estimated from laboratory experiments and micromechanical models. We use a porosity-permeability exponent n = 2, appropriate for small porosity as found in figure 7(c), from the micromechanical model of . At larger porosity, the exponent may be higher according to laboratory experiments (e.g. Wark & Watson 1998;Connolly et al. 2009;Miller et al. 2014). The sensitivity of shear viscosity to porosity was estimated experimentally to be λ * = 26 by Mei et al. (2002). Micromechanical models (Takei & Holtzman 2009a; indicate that this factor decreases with increasing porosity, but is similar to the estimate λ * = 26 when porosity is fairly small (see, e.g., Fig. 13 of Rudge (2018)). The bulkto-shear viscosity ratio has been extensively debated and has not been directly measured experimentally. Micromechanical models offer different predictions depending on assumptions about the microphysics of creep. There are two main categories. Models that assume viscous deformation at the microscale have bulk viscosity proportional to η/φ, so the bulk-to-shear viscosity ratio is very large (O(10 2 ), e.g., Simpson et al. 2010). Models that assume diffusion at the microscale (either volumetric or grain boundary diffusion) have a much weaker sensitivity to porosity and the bulk-to-shear viscosity ratio is moderate (O(1), e.g., Takei 1998; Rudge 2018). We take the estimate ζ0/η0 = 5/3 as our reference case and consider the possibility that the ratio is much higher in section 5.4.
The melt velocity depends on the maximum degree of melting, for which we take Fmax = 0.2 as a typical value. It also depends on the ratio Q0/U0. As a reference value we take Q0/U0 = 6.3×10 4 , which for U0 = 3 cm/yr corresponds to a maximum melt velocity of 4 m/yr. Q0 is sensitive to the reference permeability of the mantle, which is poorly constrained. Rees Jones & Rudge (2020) argue that the maximum melt velocity is faster than this. Here we make a relatively conservative choice such that in our reference case, both reaction and shear have the potential to make similar contributions to the growth of porosity bands. We consider faster and slower melt segregation in section 5.4. The overall amount of reactive melts generated depends on the parameter group βH/α ≈ 0.2 (probably within a range 0.15-0.3), based on previous studies of the reaction-infiltration instability (Aharonov et al. 1995;Rees Jones & Katz 2018)

Combined instability at MORs
We now use reference parameters to combine this MOR background state with the calculation of the linear growth rate. We estimate the growth of instabilities caused by reaction and shear. Figure 6(b) illustrates our approach, which is based on evolving a local parcel of porosity bands along each streamline of the solid mantle flow, as we now describe in detail.

The role of porosity advection
Our previous calculations in section 3 concerned the growth rate of instability in an infinite, uniform medium. For a midocean ridge however, we must consider how to treat the advection of porosity in equation (5). We write that equation as where f is a general source term representing the effect of compaction and reaction. This equation contains the only time derivative and the only solid-velocity advection term in the overall set of equations and its treatment has been considered by previous studies. Spiegelman (2003) showed that the rotational part of a linear flow causes an evolution in the angle of porosity bands. For simple shear, vs ∝ [z, 0, 0], Spiegelman (2003) used a generalized linear analysis where the wavevector depends on time such that φ ∝ exp(ik(t)·x). Butler (2010) showed that both a pure and also a simple shear flow cause the wavelength of porosity bands to increase. Gebhardt & Butler (2016) extended the methodology of Spiegelman (2003) to a general flow with translation and shear, and applied it to the solid velocity field of a midocean ridge setting.
In this section, we use the same approach as Gebhardt & Butler (2016). The only new aspect of our calculation involves the growth rate σ, where we account for growth of instabilities by reactive infiltration as well as shear. We now give a slightly expanded justification for the methodology. We introduce a local co-ordinate system about some arbitrary point x0 and let x be the (small) displacement from that point. We then Taylor expand the solid velocity to first order about the point using the velocity gradient tensor vs(x0 + x) ≈ v0 + ∇vs · x, v0 = vs(x0), where the velocity gradient tensor ∇vs is evaluated at x0. Thus the velocity gradient tensor in the approximation is locally a constant, so the corresponding term ∇vs · x is a linear shear flow of the type discussed in section 2. We express the porosity as a generalized normal mode: where φ0 is the background state and k(t) is the wavevector of a disturbance and s(t) determines its amplitude. The wavevector and amplitude evolve according to where D Dt ≡ ∂ ∂t +v0·∇ is the Lagrangian derivative (or equivalently the derivative along a streamline of the solid flow). These equations are sufficient to ensure that the porosity advection equation (80) is satisfied, provided f ≈ f0+σ(φ−φ0), where f0 is the background part of the source term. This can be verified by substituting equation (82) into equation (80) with the velocity expanded using equation (81). The terms involving the uniform translation v0 are accounted for by the switch to the Lagrangian derivative. The terms involving the wavevector can be shown to cancel by taking the scalar (dot) product of equation (83) with the vector x.
This approach to evolution of porosity perturbations is local, and is only valid when the solid velocity varies on some scale very much larger than the wavelength of the perturbation. It can be made precise in certain situations, including that of Spiegelman (2003), which considers a uniform background state and a linear shear flow. However, as Gebhardt & Butler (2016) discuss, for mid-ocean ridges it should be thought of as a reasonable, if ad hoc, estimate of the behaviour of porosity bands. This is because the calculation of σ was for an infinite domain with a uniform background porosity, uniform magma flow field and linear mantle shear flow. These will all be reasonable approximations provided the wavelength of the porosity bands is small, in which case the bands vary on scale much shorter than that over which the background porosity evolves, and so the growth rate calculation should be reasonably accurate.

Integration along streamlines
We now integrate equations (83,84) along streamlines of the solid flow, as shown in figures 6(b) and 7(d). It is possible to do this in Cartesian co-ordinates, but preferrable to work in polar co-ordinates because the equation of a streamline is particularly simple. Consider a streamline (r(t), ϑ(t)) that starts at (r0, ϑ0) when t = 0. By choosingr0 = 1/ cos ϑ0, all the streamlines start at the base of the melting region. Then a streamline is defined bỹ rΘ(ϑ) = constant ⇒r(ϑ) = Θ(ϑ0) cos ϑ0
Thus given ϑ, we have an explicit expression forr and hence (x,z) along the streamline, if required. Then we have an evolution equation for ϑ, namely Since we are integrating along streamlines, the material derivatives D/Dt can be replaced by full derivatives d/dt. Finally, we can change variables to make ϑ the independent variable instead of t using equation (86). Thus This system of equations can be integrated over the range of ϑ0 ≤ ϑ ≤ ϑ1 using any standard ODE solver (we use the MATLAB routine ODE45). This procedure is repeated for a range of initial positions specified by ϑ0. The system is linear in the initial amplitude so we take s = 0 as an initial condition. We also need to specify an initial condition on the wavevector (see section 4.2.4). Note that the change of variables to ϑ breaks down exactly beneath the ridge axis (by construction, since ϑ is not varying along that streamline). For that special case, it is easiest to use z as the independent variable.

Overall growth rate
The amplitude of the instability evolves along a streamline according to equation (88). We neglect the imaginary part of the growth rate (which gives rise to transient waves) and use the simplest expression for the real part of the growth rate (52), in which the contributions from reaction and shear can be computed separately. Thus we split the amplitude s from equation (82) into a part arising from reaction sreaction and a part arising from shear s shear . For reaction-driven instabilities, we find Thus the magnitude of the reactive growth rate depends on the following dimensionless parameters: the permeability exponent n = 2, the overall amount of reactive melts generated βH/α = 0.2, and the relative magma velocity w0/U0 given by equation (79). This latter ratio can be significantly greater than 1, allowing the contribution from the reactive instability to be significant. For shear-driven instabilities, we find in which we used equation (70). Crucially, this is independent of the plate half-spreading rate. Although a faster spreading rate increases the strain rate, which increases σ shear , it also increases the solid mantle velocity, thereby reducing the time spent to move through the partially molten region (proportional to H/U0). Instead, the the shear-driven instability depends only on a combination of rheological properties 2λ * / 4 3 + ζ 0 η 0 , relative position in space (via ϑ and ϑ1), and orientation of the wavevector.

Initial conditions and evolution of the wavevector
We specify the wavevector at the start of a streamline as an initial condition. The wavevector then evolves according to equation (87), where∇ũ, obtained from the corner flow, is given by equation (68). We choose an initial condition ky = 0. Then by equation (87), ky remains zero. This choice of initial condition is motivated by the fact that a mid-ocean ridge has extension in the horizontal, which favours a wavevector orientation with ky = 0, as discussed in section 3.2.
We take two approaches to specify the initial condition. First, we prescribe a single initial orientation of the wavevector. Second, we prescribe a uniform distribution of wavevector orientation, an approach also taken by Gebhardt & Butler (2016). Figure 8 shows the effect of the corner flow on the wavevector, starting from a single initial orientation. The flow acts to rotate the wavevector in a clockwise sense to the right of the ridge axis (decreasing θ). The wavenumber variation depends on whether the flow resolved along it is extensional or contractional. However, the maximum change of wavenumber (and hence of wavelength) is a factor of about 2 of its initial value. Given that the growth rate is only weakly sensitive to wavelength (assuming it is smaller than the compaction length), we can infer that the effect of solid flow on the wavenumber is mainly through rotation. Figure 8. Effect of shear on orientation and magnitude of wavevector. In each panel, the initial magnitude of the wavevector is 1 and the initial orientation is: (a) θ 0 = π/2, (b) θ 0 = π/4, (c) θ 0 = 0, (d) θ 0 = −π/4. Black line segments show the orientation of the wavevector and red line segments show the corresponding porosity bands, which are perpendicular to the wavevector. The colourscale shows the wavenumber, with red colours corresponding to increased wavenumber (reduced wavelength) and blue colours corresponding to reduced wavenumber (increased wavelength).

RESULTS: GROWTH OF INSTABILITIES AT MORs
Results are described in the following order. First, in section 5.1, we calculate the maximum possible growth of the instability by assuming that the wavevector is always instantaneously in the optimal orientation. Second, in section 5.2, we calculate the growth for a prescribed initial wavevector orientation that evolves along streamlines. Third, in section 5.3, we calculate the growth for a distribution of initial wavevector orientations. Finally, in section 5.4, we discuss how the results depend on the choice of parameters.

Orientation-independent, maximal growth
To facilitate an understanding of instability growth in the context of our mid-ocean ridge background state, we first present results where growth rates are maximised over all possible perturbation orientations θ. In particular, for reaction we take cos 2 θ = 1 in equation (89), and for shear we take G(θ, ψ) = 1 in equation (90). This approach gives an upper bound on growth, because the factors involving the wavevector orientation are always less than or equal to 1 in magnitude. Figure 7(d) compares the ratio of local maximum growth rate due to shear versus that due to reaction. Reaction is dominant beneath the ridge axis and shear is dominant off-axis along the base of the lithosphere. This reflects the pattern of background strain rate (which is zero on the axis) and magma segregation speed (which decreases away from the axis), as shown in other panels of figure 7. Figure 9 shows the maximum possible accumulated growth of perturbations. The maximum possible growth is comparable between the two mechanisms, but this is sensitive to the choice of parameters. If the magma velocity were assumed faster or the ratio of bulk to shear viscosity were higher, the reactive growth would be much larger than the shear-driven growth (see section 5.4 for a fuller sensitivity analysis). Even the reference parameter choices (that allow both reaction and shear to contribute) emphasize an important result: beneath the ridge axis, reaction dominates in the formation of channels. The dashed magenta contours in figure 9 highlight where the amplitude factor s = 7, which corresponds to a growth in amplitude of exp(s) ≈ 10 3 according to the definition in equation (82). This suggests that deep channels were probably formed through reactive, rather than shear-driven instability. This would only be reinforced by inclusion of volatile-driven reactive melting (Keller & Katz 2016).

Orientation dependence of growth
We now turn our attention to the consequences of wavevector orientation. Figure 10 shows a set of four models, plotted in terms of accumulated growth of the perturbation. In each case, we set the perturbation wavevector to have orientation θ0 as it enters the melting region from below. The wavevector then evolves along streamlines of the corner flow. The left column shows the total growth s total , which is the sum of the growth from reaction sreaction (middle column) and the growth from shear s shear (right column). Overall it is evident that accumulated growth is sensitive to the initial orientation.
Row (a) shows that when the initial wavevector is vertical (θ0 = π/2; the porosity bands are initially horizontal), the instability is suppressed. There is contraction across the bands, so the shear-driven instability has a negative growth rate. The bands remain close enough to horizontal that the reactive mode of instability only grows slowly. The net effect is that the overall growth is negative. Row (b) shows that when the initial wavevector has an angle θ0 = π/4, there is some growth of the instability. Again, there is contraction across the bands, so sheardriven instability has a negative growth rate. The bands also have a larger component in the vertical direction, so reactively-driven instability is now more significant. In total, the reaction-driven instability is large enough to offset the shear-driven suppression of the instability and give net positive growth.
Row (c) shows that when the initial wavevector is horizontal (θ0 = 0) and the porosity bands are initially vertical, the instability grows rapidly. These bands start in the orientation most favourable to reactively driven instability. However, they are rotated into an orientation with relatively large extension across bands, so the shear-driven instability is also significant. Rotation along streamlines slightly reduces the reactive growth, but this reduction is insignificant. In total, reaction and shear cooperate to drive strong growth of porosity bands. Row (d) shows that when the initial wavevector has an angle θ0 = −π/4, the instability is partially suppressed. There is some initial growth of both reactive and sheardriven instability. However, following the streamlines, the wavevector is rotated into an orientation where there is very little reactive growth and there is decay caused by contraction across the bands. Thus, in total, by the time the lithosphere is reached, the porosity bands are suppressed. Figure 11 considers the evolution of a spectrum of initial wavevector orientations. This is evaluated along two particular streamlines, one near the ridge axis and one further off-axis, shown by the dashed curves in figure 7(d). At the base of the melting region, independent of position ϑ0, we assume that the initial magnitude of the wavevector is uniformly distributed, i.e., independent of θ0. The spectrum is then evolved along streamlines. We plot the final accumulated growth and final wavevector angle for two particular streamlines, attained when these streamlines terminate at the lithosphere. Row (a) shows results for a streamline near the ridge axis. Here, the reactive growth favours wavevectors orientated very close to horizontal (vertical porosity bands). The shear-driven instability favours wavevectors orientated with slightly positive initial angle for reasons discussed above. The overall growth favours a wavevector orientation that is intermediate between these angles. The final orientation θ1 associated with the greatest accumulated growth is rotated by the shear flow. So high-porosity bands might be expected to correspond to a final wavevector orientation in the pale yellow shaded region range, roughly −π/2 < θ1 < 0.

Growth with a spectrum of initial wavevector orientations
Row (b) shows results for a streamline further from the ridge axis. The results are similar to row (a). The only important difference is that reactive growth favours wavevectors orientated with a slightly positive initial angle. This is because such bands accumulate more reactive growth as they are rotated clockwise into the vertical orientation by the shear. For this example, the overall growth favours an initial wavevector angle that is very similar to that favoured by both reaction-driven growth and shear-driven growth separately, with the peak s occurring at a very similar angle for each mechanism. Further calculations (not shown) find similar patterns of behaviour even further off-axis, so the results in row (b) are illustrative of the general pattern away from the immediate vicinity of the ridge axis (row a). Figure 12 shows the evolution of perturbations that have maximal ultimate growth over all initial orientations θ0. This maximum is evaluated for each streamline separately, with the aim of highlighting the perturbations that would be most likely expressed in a full, non-linear solution (albeit the full non-linear solution may behave differently, as discussed in section 6.2). An alternative approach is to optimize over the initial wavevector orientation independently at each point in the interior of the melting region, rather than just at the end of each streamline. We present calculations using this alternative approach in appendix B. Figure 12 demonstrates that the most unstable orientation varies only slightly with distance off the axis. At every location, the most unstable wavevector is close to horizon- Figure 10. The growth of instabilities initialized with a single wavevector orientated with row: (a) θ 0 = π/2; (b) θ 0 = π/4; (c) θ 0 = 0; (d) θ 0 = −π/4. The left column shows the total growth s total which is the sum of the growth from reaction s reaction (middle column) and the growth from shear s shear (right column). The colour-scale is the same for all panels (and the range is clipped). Black line segments show the orientation of the wavevector and red line segments show the corresponding porosity bands, which are perpendicular to the wavevector.
tal and the corresponding porosity bands are close to vertical. The calculations reinforce the point that was discussed when considering the upper bound on growth shown in figure 9 -that all channels are initially formed by the reactive mode of instability, and that axial channels are dominated by reactive growth. In this case, shear actually reduces the growth for the deeper part of the streamline. For off-axis streamlines, shear contributes to the growth of instability as the streamlines approach the lithosphere. However, the roots of channels are always associated with reactive rather than shear-driven instability. Figure 13 shows that this emphasis on the importance of reactive instability is robust to changes in the parameters considered. In row (a), the ratio of bulk to shear viscosity is increased to ζ0/η0 = 10, which is six times higher than Figure 11. The effect of the orientation of the wavevector θ 0 on the total accumulated growth s at the end of a streamline (left column) and orientation of the final wavevector θ 1 (right column). A black cross marks the most unstable initial wavevector orientation. The pale yellow shaded region highlights wavevector orientations with a total growth within 80% of the maximum. The top row (a) is an example for a streamline near the ridge axisx 0 = 0.01 and the bottom row (b) is an example for a streamline further from the ridge axisx = 0.2. These streamlines are shown as dashed magenta curves in figure 7(d).

Parametric sensitivity
the reference case. This has the effect of suppressing the shear-driven mode of instability, such that the total accumulated growth is dominated by reaction throughout the melting region. Consequently, the most unstable orientation of the wavevector is close to horizontal, since this produces the vertical, high porosity channels favoured by reaction. The solid flow still plays a role in rotating the wavevector. In row (b), we show that a similar pattern is obtained by increasing the relative magma flow speed Q0/U0 by a factor of 10, which corresponds to a maximum melt speed of 13 m/yr. This increases the reactive growth rate rather than decreasing the shear-driven growth rate, but the relative effect is the same as for variations in the viscosity ratio (note the different colour scale). The only situation that allows shear to play a greater role is when the bulk-to-shear viscosity ratio is low (as in the reference case) and the melt velocity is also low. For this case, taking Q0/U0 as 16 times smaller than the reference case, which corresponds to a maximum melt speed of 1 m/yr, row (c) shows that in most of the domain, shear is more important than reaction. Very close to the ridge axis, reaction remains dominant. However, as we discussed in section 4.1.3, this segregation rate is probably too slow to satisfy observational constraints. Thus, under more realistic choices of parameters, our model predicts that reaction-driven instability plays the dominant role.
In figure 14 we consider a reduced lithospheric slope, such as might correspond to a faster spreading rate. The results are more complex than in the reference geometry, so we plot them for a set of different initial wavevector orientation to explore this complexity (as in fig. 10). The pattern of reaction-driven growth is similar to calculations with a more steeply dipping lithospheric base. This is as expected, given that it is not directly sensitive to the strain-rate field. However, the pattern of shear-driven growth differs, especially in rows (b) and (c). In row (b), we find a positive contribution to the instability from shear as the streamlines approach the lithosphere, whereas the equivalent contribution in figure 10(b) was slightly negative. Conversely, in row (c), we find a slightly negative contribution whereas the equivalent contribution in figure 10(c) was positive. These differences reflect the different orientation of strain rate that results when the base of the lithosphere has a shallower dip. Nevertheless, the overall picture of a reactive instability that dominates the dynamics and favours sub-vertical high-porosity channels persists under these conditions.

Summary of results
The results above investigate channelized melt extraction from the mantle and the combined role of two known mechanisms of flow localisation, reaction-and shear-driven instability. The theoretical framework developed in section 2 allows us to simultaneously describe these two mechanisms in a manner consistent with published results obtained separately for the reaction-infiltration instability and the sheardriven instability. In section 3, we showed that the relative importance of shear-driven versus reaction-driven instability is governed by the dimensionless ratio given in equation (51), Figure 12. The evolution of the most unstable initial orientation of the wavevector from the bottom row of Figure 11. Other figure details are as in Figure 10.
which we rewrite here: The ratio S represents the ratio of growth rates due to shear and due to reaction. S is controlled by a particular combination of mechanical material properties the reactivity of the system β/α (which has units of m −1 ), the background rate of melt flow w0, and the background solid strain rateγ0. The dimensionless parameter S is distinct from both the "Fiji" number Φg = ∆ρg and the Damköhler number discussed in the review of Kohlstedt & Holtzman (2009). The parameter Φg is related to the ratio of shear-driven melt velocity to buoyancy-driven melt velocity and would come into the imaginary part of the growth rate. The imaginary part is associated with porosity waves; here we chose to focus on the real part of the growth rate, which controls the amplitude of porosity bands. The Damköhler number controls the degree of disequilibrium. We showed that this plays only a modest role in affecting the overall growth rate; instead, the crucial parameter is the reactivity β/α, which appears in S. The definition of S highlights the crucial role played by the material properties of partially molten rocks, some of which are not well constrained. Indeed, there remain important differences between micromechanical models. Perhaps the greatest uncertainty is the case of the bulk viscosity, as demonstrated by the contrast between model predictions (Takei & Holtzman 2009a,b;Simpson et al. 2010;. As a consequence of these uncertainties, robust, leading-order features of laboratory experiments (e.g., Holtzman & Kohlstedt 2007) such as the orientation of highporosity bands, their size, spacing and rate of emergence remain challenging to predict quantitatively (Alisic et al. 2016). Because of this gap in our knowledge, we face significant uncertainty in extrapolating between the laboratory scale and the mantle scale. The laboratory experiments are performed in closed capsules with a fixed melt fraction deforming at very high strain rates, whereas the MOR system is open to melt flow and deforms a million times more slowly.
In this context, we have opted for the simplest form of model that captures the shear-driven formation of porosity bands. In particular, we have not considered a power-law shear viscosity (Katz et al. 2006), anisotropic viscosity Qi et al. 2015), or the stabilizing influence of surface tension (Parsons et al. 2008;King et al. 2011a;Bercovici & Rudge 2016). Our methodology could be extended to include these effects, but such extensions are most worthwhile once there is a more settled and complete understanding of the shear-driven formation of porosity bands in isolation.
Our results also clarify the geometric controls on the combined instability. The reaction-infiltration instability has only one preferred direction in the context of our modelthe vertical. This is because the background magma flow direction and the solubility gradient are aligned with gravity. In the absence of shear, there is rotational symmetry about the vertical, so the instability leads to the formation of tube-shaped regions of elevated porosity. The presence of a large-scale, solid shear flow breaks this symmetry. Provided the deviatoric stress in the horizontal direction is extensional rather than contractional, even a small amount of shear favours the formation of tabular, high-porosity bands (section 3.2). This situation applies at a mid-ocean ridge, for example. The orientation of the resulting bands is intermediate between that favoured by reaction (vertical) and that favoured by shear, and is controlled by the parameter S. Thus the tabular geometry of dunite bodies in ophiolites is consistent with a reactive origin combined with extension in the horizontal direction. The favoured orientation is only weakly affected by considering the full pressure-dependence of the solubility gradient that drives the reaction infiltration instability (see appendix A, figure A1). Nonetheless, the pressure dependence of the solubility gradient may be important in the nonlinear development of channels, particularly in the presence of strong lateral pressure gradients. This suggests an important role for numerical models, as discussed in section 6.3.
We briefly discussed the dependence of the growth rate on wavelength, which is potentially significant in setting the length scale of dunites. Growth is suppressed at length scales that exceed the compaction length, so channels have a smaller scale in the direction of the wavevector than the compaction length, which is thought to be about a kilometre in the mantle. In previous work, we showed that the reaction-infiltration instability can lead to channels growing Figure 13. Sensitivity experiments showing the results with the optimal initial wavevector orientation. Row (a) has a higher bulk-to-shear viscosity ratio ζ 0 /η 0 = 10, which is 6 times larger than the reference case. Row (b) has a higher melt velocity ratio Q 0 /U 0 = 6.3 × 10 6 , which is 10 times greater than the reference case. Row (c) has a lower melt velocity ratio Q 0 /U 0 = 4 × 10 3 , which is 16 times smaller than the reference case. Note the different colour scales. Other figure details are as in Figure 12. with a scale consistent with geological observations, within the considerable parametric uncertainty (Braun & Kelemen 2002;Rees Jones & Katz 2018). Given this parametric uncertainty and the incomplete understanding of the length scales of shear-driven porosity bands, it is premature to draw definitive conclusions. We applied our model of the combined instability to a mid-ocean ridge using the method proposed by Gebhardt & Butler (2016), as described in section 4. Our results (section 5) are summarised in figure 15, where panels a-c show the accumulated growth s as a function of depth for three different corner-flow streamlines that ascend from the bottom of the melting region. The contribution of shear (green) is significant only at shallow depths and/or far from the ridge axis. Geochemical evidence requires channels to form at least 15 km beneath the Moho (Kelemen et al. 1997) and U-series disequilibrium and reactive-flow models suggest it is possibly much deeper Keller & Katz 2016;Liu & Liang 2019). Together this indicates that instantaneous growth rates are dominated by reaction and affected by shear only along the base of the lithosphere.
Nonetheless, the shear associated with tectonic-scale flow is important in that it promotes the tabular geometry of channels and affects the orientation of channels by rotating their tops away from the ridge axis, which corresponds to a decrease in wavevector angle along the streamline, as shown in figure 15d. Most of this rotation happens very close to the upper end of a streamline, where the streamline turns the corner and terminates at the base of the lithosphere. The most unstable initial orientation of bands is close to vertical, with their tops tilted slightly towards the ridge axis. Bands with this orientation grow rapidly, mainly due to the reactive instability. Near the end of the streamline, the bands are rotated further and can end up lying between the vertical and the horizontal. We explored whether these results are robust to the choice of parameters used, including the angle from the horizontal to the base of the lithosphere. The only combination of parameters that favoured shear over reaction is a very slow melt velocity ( 1 m/yr) with a small bulkto-shear viscosity ratio (O(1)). There is evidence for much more rapid melt extraction ( 10 m/yr) based on the Icelandic deglaciation (e.g. Jull & McKenzie 1996;Maclennan et al. 2002;Eksinchol et al. 2019;Rees Jones & Rudge 2020) and U-series data (e.g. Iwamori 1994;Kelemen et al. 1997;Jull et al. 2002;Stracke et al. 2006;Elliott & Spiegelman Figure 14. Sensitivity experiment with a shallower dip π/12 of the lithospheric base showing a range of initial wavevector orientations. Other figure details are as in Figure 10. 2014). Given much more rapid melt extraction, it is most likely that reaction is dominant over shear in general.

Limitations of analysis
Our calculations, like those of Gebhardt & Butler (2016), rely on a separation of scales between the sub-compactionlength scale of porosity localisation ( 1 km) and the tectonic scale of the mid-ocean ridge (∼ 100 km). This separation arguably enables us to embed our idealised calculation of the local growth rate from section 3 in tectonic models of mid-ocean ridge magmatism. However, we have not shown that these scales are truly separate and, indeed, there are reasons to question this. For example, local perturbation growth could create a pattern of anisotropic material properties (e.g., permeability, viscosity) that would feedback on the large-scale dynamics. It is also worth noting that the same scale separation presents a severe challenge for numerical models that discretise a two-or three-dimensional space. To capture the interaction between the large and small scales requires either extremely high grid resolution over a very large domain or a more sophisticated, multi-scale approach (e.g., Kevrekidis et al. 2003).
The linearity of the current approach is another important limitation. The linearized equations apply rigorously only to the initial stages of channel formation; their formal validity breaks down as the perturbation to the background state grows. Exponential growth cannot continue indefinitely and, instead, the channel amplitude saturates (Spiegelman & Kelemen 2003;Liang et al. 2010;King et al. 2010). At the same time, the lithological imprint of channelised flow (replacive dunites) is advected and rotated by the mantle flow. Lithological structure may serve to lock in the pattern of reactive channels and limit overprinting by shear-driven growth in contrasting directions. This indicates that two-dimensional numerical models have an important role to play in understanding the nonlinear evolution of channels, channel coalescence and the overall arrangement of magmatic localisation. Unfortunately, such numerical models sometimes fail as the localisation becomes more pronounced (e.g. Vestrum & Butler 2020), perhaps reflecting our incomplete understanding of the physics of partially molten rocks. Numerical diffusion may also mask behaviour that is expected based on the results obtained here. In one possible example of this, Katz (2010) found that shear-driven porosity bands did not emerge in two-dimensional simulations of melt transport beneath a mid-ocean ridge, even with exaggerated viscous weakening by porosity. However, Katz & Weatherley (2012) demonstrated that inherited lithological heterogeneity in the mantle can lead to sharp localisation of melting and melt transport. And while the chemical heterogeneity imposed in that case may be extreme, models by Keller et al. (2017) predicted the emergence of channelized flow due to volatilerich flux melting beneath a mid-ocean ridge. Our approach based on linearization of the governing equations complements these numerical studies.
Other potentially significant issues relate to the background state about which we linearize. It is probably a good approximation to consider that the solid flow is only minimally affected by the magma flow provided the porosity remains small. However, we used a simple Newtonian viscosity for the solid flow; this could be extended to consider a non-Newtonian or anisotropic viscosity. These effects alter the behaviour of the shear-driven instability, reducing the angle of porosity bands. Nevertheless, the principle remains that the bands grow most rapidly at an angle between that favoured by reaction and that favoured by shear. We also assumed that the magma flow was purely vertical, driven by buoyancy. While this can be the dominant contribution to magma transport, there are other mechanisms (discussed in the following section) that give rise to a more complex, not purely vertical, pattern of magma flow at mid-ocean ridges.
Another related issue that could modify our predictions for mid-ocean ridges is the background compaction rate associated with melt segregation. If the compaction viscosity is a decreasing function of porosity, background compaction acts to stabilize the system against exponential perturbation growth. Hewitt (2010) showed, in the context of a meltingcolumn model, that this effect is significant if the bulk viscosity has an inverse dependence on porosity. The stabilising effect is much weaker if the compaction viscosity varies only logarithmically with porosity Rees Jones & Katz 2018). This could be reassessed in the mid-ocean ridge geometry.

Implications
A leading-order observation about mid-ocean ridges is that the volcanic zone is much narrower (∼10 km) than the lateral extent of the partially molten region (∼100 km), which means that melt must be focused laterally towards the ridge axis. Spiegelman & McKenzie (1987) and Morgan (1987) hypothesised that dynamic pressure gradients suck melt towards the ridge. Sparks & Parmentier (1991) and Spiegelman (1993) proposed that melts migrate to the ridge through a sub-lithospheric decompaction channel. More recently, Turner et al. (2017) and Sim et al. (2020) have argued for 'melting-pressure focusing' associated with gradients in compaction pressure. In addition to these mechanisms, Katz et al. (2006) suggested that shear-driven porosity bands create an effectively anisotropic permeability, and that they have an orientation such this anisotropy focuses melt toward the ridge (see also Morgan 1987;Daines & Kohlstedt 1997;Kohlstedt & Holtzman 2009;Liu & Liang 2019). In contrast, the present calculations indicate that high-porosity channels typically have a sub-vertical orientation in the region moderately close to the ridge (see figure 13 and 15d). Beneath the lithosphere, they remain sub-vertical but tend to rotate away from the ridge axis due to the corner flow. The anisotropic permeability structure that might arise from this pattern would not contribute much to melt focusing. On the other hand, melt focusing by pressure gradients might lead to reaction-induced channels that point toward the ridge axis, aligned with the magmatic flow direction (Rabinowicz & Ceuleneer 2005). Numerical models could be used to test the hypothesis that lateral pressure gradients additionally drive the formation of diagonal reactive dissolution channels through pressure-dependence of the solubility gradient (appendix A), thereby enhancing focusing (Spiegelman et al. 2001).
Coherent alignment of melt within the mantle could give rise to anisotropy of seismic wavespeeds. Measured anisotropy might therefore be related to the predictions above if the influence of melt can be disentangled from other causes of anisotropy. Kendall (1994) and Blackman & Kendall (1997) recognised the potential for grain-scale alignment of melt to shape the pattern of anisotropy beneath mid-ocean ridges. Later, Holtzman & Kendall (2010) argued that localised melt-fraction perturbations (channels or bands) would also induce seismic anisotropy. To explain observations of seismic anisotropy beneath some mid-ocean ridges, Holtzman & Kendall (2010) and Nowacki et al. (2012) invoke sheets of higher melt fraction sub-parallel to the lithosphere-asthenosphere boundary (LAB) at some 50 km from the ridge axis to create a tilted transverse isotropy (TTI). While such alignment was predicted by Katz et al. (2006) and Holtzman & Kendall (2010) on the basis of the shear-driven instability, results presented here cast doubt on it and hence on the TTI hypothesis. However, plate spreading at the Main Ethiopian Rift is associated with a fast seismic direction that is ridge-parallel (Kendall et al. 2005;Hammond et al. 2014). The sub-vertical, tabular magmatic structures that we predict to arise from the shear-modified reaction-infiltration instability are consistent with such anisotropy. Alignment of magmatic features subparallel to the shallowly dipping LAB may instead arise by the dynamic response to magmatic flow toward an imperme-able boundary (Sparks & Parmentier 1991;Hewitt & Fowler 2008). Alternatively, seismic anisotropy may arise from a grain-scale texture (e.g., lattice-or melt-preferred orientation (Holtzman et al. 2003b;Qi et al. 2018)) rather than accumulated growth of macroscopic localisation patterns, as initially proposed by Kendall (1994).

Conclusions
• We developed a consistent framework to model the combined growth of reaction-driven and shear-driven instabilities in the partially molten mantle and calculated their linear growth rate. We applied that framework to the melting region beneath mid-ocean ridges.
• The reactive-infiltration instability is dominant over most of the melting region and, in particular, along its base and close to the ridge axis, where mantle flow is vertical. This gives rise to sub-vertical, high-porosity channels that form sufficiently deep within the melting region to explain the observed disequilibrium between erupted lavas and the harzburgitic upper mantle.
• The presence of even a small amount of horizontal extension favours tabular channels over tube-shaped channels, consistent with the morphology of observed dunites.
• The shear-driven instability may contribute to further growth of channels along the base of the lithosphere. However, the orientation of these channels is set by prior growth in the reactive-flow regime and hence they remain subvertical, despite some rotation by the corner flow.
• Within the limitations of our study, shear-driven melt bands do not function as a mechanism for melt focusing at mid-ocean ridges. Moreover, we predict that bands are not aligned consistently with the shallow-dipping seismic anisotropy that has been obtained by inversions.

ACKNOWLEDGMENTS
Insightful comments by P. Kelemen and an anonymous reviewer helped us to improve the manuscript. D.W.R.J. acknowledges research funding through the NERC Consortium grant NE/M000427/1. The research of R.F.K. leading to these results received funding under the European Union's Horizon 2020 research and innovation programme, grant agreement number 772255. We thank the Deep Carbon Observatory of the Alfred P. Sloan Foundation.

Data availability statement
Software that implements the theoretical methods presented this article is archived by Zenodo and available at https://doi.org/10.5281/zenodo.4618182 (Rees Jones et al. 2021).

APPENDIX A: PRESSURE-DEPENDANT REACTIVE MELTING
The reaction-infiltration instability is driven by a chemical solubility gradient. In the main text, we assumed that this gradient was vertical. However, chemical solubility depends not on depth directly but rather on pressure. If the pressure is dominantly lithostatic, then the solubility gradient will be vertical, motivating our assumption in the main text. The purpose of this appendix is to investigate the potential role of pressure-dependent melting by relaxing the assumption that the thermodynamic pressure is dominantly lithostatic. The shear-driven instability creates a pressure gradient that can, in turn, feedback on the reactive instability through the pressure-dependent reactive melt rate. So this appendix allows us to consider another potential mode of coupling between the two types of instability.
We redo the linear stability from section 2 to account for the full pressure-dependence of the reactive melt rate. To focus on this effect, we strip out other parts of our earlier analysis from the outset. These assumptions led to the simplified growth rate reported in section 3.1, so we compare with results reported in that section.

A1 Models for pressure-dependent reactive melt rate
We start with the model for the reactive melt rate originally given in equation (14), which we replace by vD · ∇ceq = αΓ, (A.1) where we repeated our earlier simplification (φv l ≈ vD).
The solubility gradient ∇ceq plays the role of ∇(βz) earlier. We next assume that the solubility gradient is linearly related to the pressure gradient. Over the full-depth of the melting region this is a simplification, but will be valid locally, consistent with the approach taken in section 2. Thus where m > 0 is a proportionality constant and ∇P therm denotes the thermodynamic pressure gradient. The concept of thermodynamic pressure needs to be carefully considered for systems out of equilibrium. Jull & McKenzie (1996) consider two main possibilities, ∇P therm = ∇P l and ∇P therm = ∇(P l − ζC), in which ∇P l is given by equation (10): Jull & McKenzie (1996) argue in favour of the second possibility, which corresponds to −1/3 times the trace of the solid stress tensor. For the purposes of this appendix, the distinction between these definitions has only a marginal effect on the results, so we consider both possibilities for completeness. Next, we define the lithostatic pressure gradient where the approximation ρ ≈ ρs is consistent with the assumption φ 1 we made in section 2. Then we define the non-lithostatic part of the liquid pressure gradient as ∇Pm = 2Ds · ∇η + η∇ 2 u + 4 3 η∇C + ∇(ζC), If the thermodynamic pressure is dominated by the lithostatic contribution, then we can combine equations (A.2) and (A.4) to obtain ∇ceq = mρsgẑ, which is consistent with the simplified approach in the main text provided Finally, we consider the solubility gradient under the two potential definitions of thermodynamic pressure above. If ∇P therm = ∇P l , then The dimensionless ratio ∇Pm/ρsg is the additional melting factor arising from consideration of the pressure-dependence of the solubility gradient. If ∇P therm = ∇(P l − ζC), then

A2 Revised linear stability analysis
We now redo the linear stability analysis with the generalized reactive melting model. We first rewrite equation (18)  We next expand out the pressure gradient term, making the same simplification that C0 = 0 as in the main text, ∇P m = ∇P l , = 2D0 · ∇η + η0∇ 2 u + 4 3 η0 + ζ0 ∇C , = ∇ψ + 4 3 η0 + ζ0 ∇C , = −2λ * η0γ0∇ ψ + 4 3 η0 + ζ0 ∇C . using the definition of δ in equation (35) and Λ in equation (36). So v D = nw0ẑφ + Λ∇ ψ − δ 2 ∇C . (A.14) As an aside, this gives another interpretation of Λ as a factor that controls the perturbed melt flow arising from the sheardriven instability. We next linearize the generalized reactive melt rate equation (A.1), Note that the base-state solubility gradient is βẑ, independent of the additional terms arising from the pressure-dependence of the solubility gradient. The perturbation to the melt rate is given by v D · βẑ + K0∆ρgẑ · ∇c eq = αΓ . (A.15) The extra terms from the pressure-dependence of the solubility gradient depend on the definition of thermodynamic pressure as discussed above. If ∇P therm = ∇P l , we linearize equation ( The dashed curves correspond to R = 1.1 and hence account for the dependence of the solubility gradient on the dynamic pressure. For the reactive instability only (S = 0, yellow curves), the preferred wavevector angle is always θ = 0 corresponding to vertical channels, independent of R.
∇P therm = ∇P l , R is only 10% different from the purely lithostatic case. Using ∇P therm = ∇(P l − ζC), reduces this already small difference further, since ζ rel < 1 and so R is closer to 1. Figure A1 shows that the angular dependence of the growth rate from equation (A.24) is only minimally affected by the choice of R. In particular, the most unstable orientation satisfies tan 2θ = 2S sin 2θe R + 2S cos 2θe , (A.26) which is very weakly affected by R. At small S, the favoured orientation is close to kz = 0 (θ = 0) so the extra stabilization due to R > 1 is unimportant, as shown by equation (A.24). At large S, the growth rate is dominated by shear, so the extra stabilizing of the reactive part of the growth rate is again unimportant, as shown by equation (A.26). Physically, the pressure-perturbations associated with the shear instability are coupled to the reactive instability through the pressure-dependence of the solubility gradient. However, this coupling only affects the imaginary part of the growth rate. The real part of the growth rate is also affected by pressure perturbations associated with the reactive instability (not the shear-driven instability). Both of these effects are relatively small, because ∆ρ/ρs 1. This analysis is only preliminary for two reasons. First, the consideration of terms ∆ρ/ρs goes beyond the Boussinesq approximation made throughout this study (see section 2.1). Second, perhaps more importantly, our study is restricted to a linear stability analysis with a vertical base upwelling of magma. So reactive melting is driven by the vertical part of the solubility gradient, which is dominantly lithostatic. In a fully developed nonlinear state, lateral pressure gradients could drive reactive melting and this might be important in understanding the coalescence of channels (Spiegelman et al. 2001).

A3 Further considerations for MORs
The importance of the non-lithostatic solubility gradient depends on its magnitude relative to the lithostatic pressure gradient. The latter scales like ρsg.
Pressure gradients arising from buoyancy-driven flow scale like ∆ρg, so their relative contribution scales like ∆ρ/ρs -the same scaling we identified in the previous section.
The viscous corner flow (section 4.1.1) gives rise to a dynamic pressure gradient. This can be estimated from Stokes equation for an incompressible fluid of constant viscosity ∇P = η∇ 2 u. (A.27) The pressure gradient is singular at the corner (x = z = 0 on the diagram in figure 6a) and decreases rapidly away from this point. To estimate the scale of this dynamic pressure gradient ∇P dyn , we scale equation (A.27) as follows. We estimate η ∼ η0, u ∼ U0, ∇ ∼ 1/L dyn , where L dyn is the distance from the corner. Then where we used the rough estimates η0 = 10 19 Pa s, U0 = 10 cm/yr, ρs = 3 × 10 3 kg/m 3 and g = 10 m/s 2 . Note that this estimate is factor of (∆ρ/ρs) 1/2 smaller than the melt focussing length scale estimated by Spiegelman & McKenzie (1987). In conclusion, within a distance of about a kilometre from the corner, the dynamic pressure gradient is larger than the lithostatic pressure gradient. However, outside this region, the lithostatic pressure gradient is much larger than the dynamic pressure gradient (the latter decreases with the inverse square of distance from the corner). Therefore, we can neglect the contribution of dynamic pressure gradients almost everywhere beneath MORs. Finally, we estimate the scale of compaction pressure gradients. Turner et al. (2017) and Sim et al. (2020) have argued for 'melting-pressure focusing' associated with gradients in compaction pressure and we base our estimates on these studies. The base state compaction rate C0 scales like the base state melting rate Γ0 (ignoring a minus sign), which we estimated in section 4.1.2. By combining estimates in that section, we find To estimate the compaction pressure gradient, we must consider the length scale over which the compaction pressure gradient varies. For a triangular melting region an appropriate scale is the depth H, since the width also scales like this.
In practice, the length scale is probably somewhat smaller than H (Sim et al. 2020 where we used the rough estimates Fmax = 0.2, ζ0 = 10 20 Pa s, and H = 60 km in addition to the estimates used in equation (A.29). Thus the compaction pressure gradient is a very small fraction of the lithostatic pressure gradient. The same conclusion can be reached by comparing figure 7 of (Sim et al. 2020) with the lithostatic pressure ρsgH ≈ 2 × 10 3 MPa. These arguments do not mean that dynamic and compaction pressure gradients can be neglected when considering melt flow (after all the lithostatic pressure gradient does not drive any flow), only that they can generally be neglected when considering thermodynamic pressure.

APPENDIX B: ALTERNATIVE APPROACH TO WAVEVECTOR OPTIMIZATION FOR MORs
For the main MOR results (section 5), we optimized the initial wavevector orientation to maximize the total growth accumulated at the end of each streamline. This is relevant if we focus at the location where geological observations of dunite channels are most readily made. However, if we focus on the fabric at depth, it might make more sense to ask what initial wavevector orientation maximizes the total growth at each point in space separately. The results are a little harder to interpret, because the corresponding initial wavevector will not be consistent along each streamline. Figure A2 shows the results of this alternative optimization procedure at three different melt segregation speeds. We also plot in the right column the optimal wavevector orientation based on equation (55), which reflects the local growth rate. The preferred wavevector orientation is always tilted in a slightly negative orientation (meaning that porosity bands are tilted slightly away from the ridge axis. Near the ridge axis (and across a wide part of the domain in the case of fast melt segregation), the wavevector angle is close to zero and porosity bands are close to vertical. This is because reaction is the dominant mode of instability here. Elsewhere, the wavevector is tilted over at angles up to about −π/4, with more negative values at slower melt segregation rates. This reflects the relative importance of shear in this case. However, for slower melt segregation, the amplitude of the porosity bands is not very high, especially at depth. The locally most unstable orientation has a similar but distinct pattern (comparing the middle and right columns). It is important, therefore, to track the accumulation of growth and the rotation of the wavevector along streamlines because the most unstable orientation at a particular point cannot be inferred purely from the local behaviour there. Figure A2. Alternative approach based on optimizing the total growth at each interior point separately, not just at the end of a streamline. The left column shows the amplitude s, the middle column shows the wavevector orientation, the right column shows the optimal orientation in terms of the local growth rate based on equation (55). We report results at a range of different melt segregation rates. Row (a) has a higher melt velocity ratio Q 0 /U 0 = 6.3 × 10 6 , which is 10 times greater than the reference case. Row (b) is the reference case. Row (c) has a lower melt velocity ratio Q 0 /U 0 = 4 × 10 3 , which is 16 times smaller than the reference case. Note the colour scale in (a) is clipped. For the middle column, we added dashed magenta contours of s total = 7. This gives an indication of the level below which the perturbations are very small.