Upscaling from mineral microstructures to tectonic macrostructures

SUMMARY

Earth’s plate tectonic behaviour arises from lithospheric ductile weakening and shear-localization. The ubiquity of mylonites at lithospheric shear zones is evidence that localization is caused by mineral grain size reduction. Most lithospheric mylonites are polymineralic, suggesting that the interaction between mineral phases by Zener pinning promotes grain size reduction and weakening. Yet this interaction only occurs where mineral phases mix at the grain scale. Phase mixing and its effect on microstructure and strength have been shown in deformation experiments and natural field samples. Our theory for the interaction between phase mixing (treated as a stress driven diffusion) with two-phase grain damage has been compared to lab experiments. But using processes at the tiny grain-scale embedded within the small hand-sample and lab scales to model large-scale lithospheric processes, requires an upscaling scheme that captures the physics from micro- to macrostructures. For example, weakening from grain-damage in zones of mixing can lead to banded viscosity structure at the small scale that manifests as viscous anisotropy at the large scale. Here we provide a new framework for self-consistently upscaling from microscopic (grain) scales, to mesoscopic (petrological heterogeneity) scales to macroscopic (tectonic) scales. The first upscaling step models phase mixing and grain size evolution in a small ‘mesoscopic’ lab-scale volume or ‘patch’, which is equivalent to a point in the macroscopic space. Within this mesoscale patch, stress driven diffusive mixing is described by an analytical solution for mineral phase fraction, provided a minimalist Fourier representation of phase fraction, and a transformation to the patch frame of reference as well as to the principal stress directions at that point. The orientation and volume fraction of mixed-phase regions can then be extracted from the analytical solution for phase fraction. The grain size and viscosity in the mixed bands are determined by two-phase grain-damage theory; the unmixed zone properties follow from mono-phase grain damage theory. The mesoscale banded viscosity field leads to a macroscale anisotropic viscosity at that point in space. But, the evolution of properties at each macroscale point involves tracking only a few quantities (phase fraction, grain sizes) rather than modelling each patch of mesoscale space as its own 2-D or 3-D system. For the final upscaling, the anisotropic viscosity field is used in a macroscale lithosphere flow model. We show an example of this scheme for a lithospheric Rayleigh–Taylor drip driven by ridge-push compressive stress, which can cause anisotropic weakening via grain mixing and damage that may help initiate subduction and passive margin collapse.

1 INTRODUCTION

Most of the prominent tectonic deformation on the solid planets and moons is manifested as localized strain in the crust and lithosphere (Thatcher 2009; Kohlstedt & Mackwell 2010; Harris & Bédard 2014; Byrne et al. 2021; Karagoz et al. 2022; Hanmer 2023). Earth’s weak plate tectonic boundaries, which span the entire breadth of the lithosphere and cut across both brittle and ductile zones, are the archetype of such deformation (for reviews and recent examples, see Bercovici 2003; Bercovici & Ricard 2014; Gueydan et al. 2014; Bercovici et al. 2015; Karato & Barbot 2018; Mulyukova & Bercovici 2019, 2023; Becker & Fuchs 2023). Although crustal and lithospheric deformation can (and often is) modelled geodynamically as a viscous or plastic (or viscoplastic) medium (e.g. Tackley 2000; Foley & Becker 2009), that treatment belies the important control of microstructural evolution—through the motion of defects like vacancies, dislocations and grain boundaries within minerals (e.g. Ranalli 1995; Karato 2008; Mulyukova & Bercovici 2018a, 2022). Microstructural response not only weakens rocks, but provides memory of deformation states (e.g. Bercovici & Ricard 2014; Fuchs & Becker 2021, 2022), can alter the transient behaviour of Earth’s materials (e.g. Montési & Hirth 2003; Montési 2004; Kidder et al. 2016; Mulyukova & Bercovici 2018a, 2022), and strongly influences mantle circulation and its correspondence to seismological observations (Behn et al. 2009; Wada et al. 2011; Dannberg et al. 2017). While there are various forms of micromechanical behaviour influencing rock deformation, mylonites especially are a tell-tale sign of the role of grain size in shear localization in the crust and lithosphere. Mylonites are metamorphic rocks that are widespread in ductile deformation zones and plate boundaries, and whose miniscule grain sizes are associated with intense deformation (e.g. White et al. 1980; Jin et al. 1998; Furusho & Kanagawa 1999; Warren & Hirth 2006; Skemer et al. 2010; Linckens et al. 2011; Hansen et al. 2013; Boneh et al. 2021). The positive feedback between deformation induced grain comminution and grain size-dependent rheology can provide the conditions for localization, whereby deformation produces weakening, which focusses deformation, which promotes more weakening and so on.

Mylonitization, especially in the lithosphere, occurs where mineral phases mix with each other (e.g. Herwegh et al. 2005, 2011; Linckens et al. 2011, 2015). Mixing of phases induces Zener pinning of one mineral phase against the other, wherein one mineral blocks the grain-boundary migration of the other mineral; this impedes grain growth and likely facilitates grain damage and comminution by dynamic recrystallization (DRX, see Bercovici & Ricard 2012; Mulyukova & Bercovici 2019). Mylonitization in monominerallic media is not nearly as pronounced, as is also the case for polyminerallic rocks comprised of unmixed phases existing as isolated monophase tectonites. The interaction between different minerals that promotes localization requires phase mixing during deformation. Phase mixing, and its effect on rock strength, have been the subject of various laboratory, field and theoretical studies (Bercovici & Ricard 2016; Tasaka et al. 2017a, b, 2020; Cross & Skemer 2017; Bercovici & Skemer 2017; Bercovici & Mulyukova 2018; Wiesman et al. 2018; Cross et al. 2020; Bercovici & Mulyukova 2021; Wiesman et al. 2023a, b; Bercovici et al. 2023; Billings & Skemer 2024). Whether ductile grain-mixing at lithospheric conditions (i.e. without cavitation or cataclasis) occurs through mechanical mixing (e.g. the ‘tooth’ model of Bercovici & Skemer 2017) or by dissolution and chemical diffusion of a metal oxide unit, like |$\rm MgO$| (Tasaka et al. 2017b), the mixing process appears to be diffusively dominated and is driven along compressive stresses and perpendicular to the interphase boundary between mineral units. These mechanisms can be treated agnostically with a generic theory for stress driven diffusive mixing (Bercovici & Mulyukova 2018, 2021; Bercovici et al. 2023).

In addition to grain size reduction, deformation of crystalline materials can change their mechanical properties through formation of anisotropic fabric. Although anisotropic fabric development in the crust has various causes, in the mantle and lithosphere, development of crystal preferred orientation (CPO) in olivine can weaken the mineral by up to an order of magnitude in the direction of net integrated strain (Durham & Goetze 1977; Zhang & Karato 1995; Kaminski & Ribe 2001; Hansen et al. 2012; Király et al. 2020). While it takes less strain to develop CPO fabric, compared to grain size reduction through DRX, the drop in viscosity in the case of the former is modest relative to that of the latter, which can be multiple orders of magnitude (Poirier 1980; Karato et al. 1980; Van der Wal et al. 1993; Bercovici et al. 2001b; Bercovici & Ricard 2012; Mulyukova & Bercovici 2017). Indeed, in the exposed lithospheric shear zones, CPO microfabric is primarily seen in the coarse-grained and less deformed parts of the rock (Ebert et al. 2007; Linckens et al. 2011), compared to the extremely deformed fine-grained regions (e.g. Warren & Hirth 2006). It is thus reasonable to assume that while CPO development influences the mechanics in the early stages of deformation, weakening through grain damage has the dominant impact at larger strains, as typical in mylonites.

Microstructurally induced weakening, be it CPO or grain damage, persists even after deformation ceases, since erasing the microfabric (e.g. randomizing crystal orientations in the case of CPO, or grain growth in the case of grain size reduction) is slow relative to the rates of tectonic plate motions. Thus, microstructurally weak portions of the lithosphere can be inherited and reactivated (Bercovici & Ricard 2014; Tommasi et al. 2009; Király et al. 2020). Regardless of the microstructural mechanism, some of the most notable effects of anisotropic viscosity on geodynamic processes include higher growth rates of gravitational instabilities, such as that involved in lithospheric delamination or collapse of passive margins (Lev & Hager 2008; Perry-Houts & Karlstrom 2019; Bercovici & Mulyukova 2021).

Deformation in polyminerallic media can also lead to viscous fabric or anisotropy in two ways, even neglecting anisotropy associated with crystal lattice orientation. First, as a medium is deformed, the heterogeneity of unmixed mineral phases (i.e. petrological heterogeneity) becomes distorted and sheared, and regardless of the original structure of the heterogeneity, the mineral phases and the boundaries between them are aligned by the shearing motion to create anisotropic fabric, depending on the strength contrast between phases (Treagus 2002, 2003; Fletcher 2004; Montési 2007; Dabrowski et al. 2012). If the weaker phase becomes interconnected and reaches a percolation limit, then the aggregate material can undergo a rapid strength loss with increasing strain (e.g. Montési 2007; Girard et al. 2016). Secondly, if phase mixing occurs during deformation as a diffusive process, its effective diffusivity is anisotropic since it is a function of the stress tensor (Bercovici & Mulyukova 2018, 2021; Bercovici et al. 2023). Thus even in an isotropic heterogeneous phase field, diffusive mixing has a preferred direction dictated by the stress, leading to mixing fronts with preferred orientation (mostly normal to compressive stress directions); the aligned mixed zones are where two-phase grain-damage and weakening occur and thus appear as bands of weakness that in aggregate provide anisotropic viscosity (Bercovici & Mulyukova 2021). These two processes can obviously interact as the shearing and alignment of weak bands will work on mixed zones as well as on the mineral phases themselves. The development of viscous anisotropy might play an important role in weakening the lithosphere to promote, for example, subduction initiation (Király et al. 2020; Bercovici & Mulyukova 2021).

Microstructural behaviour and phase mixing likely control the strength and evolution of the lithosphere, but they are all processes that occur at very different length scales. The extrapolation from microstructural physics to mantle and lithospheric physics is a long-standing problem in geophysics, and canonically treated by using a microstructurally informed, empirical and parametrized rheology in a macroscopic geodynamics model; for example the dislocation creep rheology washes away any knowledge of dislocations into an empirical pseudo-plastic, steady-state-creep law with power-law dependence of strain rate on stress. Transient behaviour or inheritance of weakness can be treated with a state variable as a proxy for microstructural defect evolution, for example, in mean grain size or dislocation density (e.g. Ricard & Bercovici 2009; Bercovici & Ricard 2012; Mulyukova & Bercovici 2022). However, while these extrapolations connect microstructure to macrostructure, none of them make the link through the mesoscopic scale of petrological heterogeneity, where phase mixing (even at its small scale) is an intrinsically multidimensional process (e.g. advection-diffusion in 2-D or 3-D), but its occurrence is crucial for promoting microstructural processes that are influenced by interphase boundary density, like grain-damage (Bercovici & Ricard 2012; Bercovici et al. 2023), as well as phase-boundary sliding (Zhao et al. 2019) and phase-boundary chemical transport and storage (Cukjati et al. 2019).

Here we propose a framework to connect these disparates scales of the lithospheric system, from microscopic (grain size), to mesoscopic (petrological heterogeneity) to macroscopic (mantle and lithosphere). Our goal is to provide rapid and efficient treatments of time-dependent evolution at the micro- and mesoscale that can be used to determine mechanical response at the macroscale, and without having to solve complicated, multidimensional and computationally prohibitive processes at all three scales.

2 UPSCALING FRAMEWORK

2.1 Overview of upscaling strategy

As already noted, we have three scales in our system, the macroscale (tectonic), the mesoscale (petrological heterogeneity) and the microscale (grain size, see for example Fig. 1). The characteristic lengths of each scale are as follows:

• The macroscopic domain is of characteristic size L, which is of order kilometres or more, although the domain may be of length L, and height H in two dimensions (and additionally a depth D in three dimensions).

• The mesoscale domain has a length scale |$\delta\ll L$|⁠, which is for example centimetres across. Each point in the macroscopic domain is a mesoscale patch of characteristic area |$\delta^2$| in 2-D, or a volume |$\delta^3$| in 3-D.

• Within the mesoscale domain are grains of characteristic size R, which are millimetres to microns in size.

Our general upscaling (and downscaling) approach is outlined below, and also graphically summarized in the flow chart of Fig. 2.

1. Begin at the macroscale to get stress and flow fields

   • We initiate the model with a simple solution for the macroscale stress, velocity and strain-rate fields (convective, or due to ridge-push or slab pull, etc.).

   • These initial stress and flow fields are typically the solution for the model with constant viscosity (i.e. before viscosity is affected by meso- and microscale processes).

2. Down to the mesoscale for phase mixing

   • The macroscopic stress and flow fields are imposed at the mesoscale to determine phase mixing

   • Grain-mixing diffusion is calculated within each mesoscale patch at each point in macro-space (i.e. one macroscale coordinate point is a mesoscale patch).

   • Within each patch we make two coordinate transformations, one to its mean velocity frame, and a second into the principal stress frame; in this frame, stresses and strain-rates are purely normal (compression and tension) and velocity corresponds to uniform pure shear and rotation.

   • The diffusive mixing model therefore has a simple analytic solution that can be calculated quickly at each mesoscale patch (and hence for every gridpoint in the macroscale domain).

   • The diffusive mixing solution is used to determine the orientation and volume fraction of mixed zones or bands, which are later used to infer anisotropic viscosity structure.

3. Down to the microscale for grain-damage:

   • The mixing bands, between monophase units, undergo two-phase grain damage with Zener pinning, which drives mylonitization and weakening in the mixed bands.

   • The monophase unmixed zones (outside the mixed bands) undergo single-phase grain damage, which leads to protomylonitzation, that is larger grains dominated by dislocation creep that are not as weak as those in the mixed bands.

   • Given the calculated grain sizes (as well as the imposed macroscopic strain-rate), we determine the viscosities in the mixed and unmixed zones.

4. Back up to the mesoscale: Viscosity structure and fabric of a mesoscale patch

   • Assuming that bands of mixed (and low viscosity) and unmixed (higher viscosity) regions form an equivalent bilaminate structure, we construct a transverse isotropic viscosity field for that patch (i.e. high viscosity for normal stresses acting on the bands, and low viscosity for shear stresses parallel to the bands) in the rotated frame of reference for that patch.

   • When transformed back to the reference coordinate frame, the transverse isotropic viscosity field leads to an anisotropic viscosity for that patch.

5. Back up to the macroscale:

   • The compilation of anisotropic viscosity values for each patch leads to a complete anisotropic viscosity field for the entire domain.

   • The anisotropic viscosity is used to calculate the new velocity and hence stress and strain-rate fields, given some forcing that drives deformation.

   • The new velocity, stress and strain-rate fields in the macroscopic scale are then used to update the mineralogical heterogeneity (phase distribution) and phase mixing in each mesoscale patch, which in turn dictates the microscale grain evolution in mixed and unmixed regions, and so on.

In the following sections we elaborate on this framework in detail.

2.2 From a given stress and flow field, to the mesoscale patch

The process of scaling up and down does not necessarily have a starting or ending point since the scheme is cyclic and iterative. But since all stories have to begin somewhere, we start with a given macroscopic flow field with velocity |${\mathbf {v}}$|⁠, associated strain-rate tensor |$\dot{\underline{\mathbf {e}}}$| and stress tensor |${\underline{{\boldsymbol{\tau }}}}$|⁠. These fields could be associated with any fluid flow and deformation in the mantle and lithosphere driven in the end by some version of convective forcing. The initial stress and flow fields can be based on simple solutions for an isoviscous system sans grain-damage and mixing at the various subscales. After we use these fields to evolve the system at the meso- and microscales, and feed that evolution back to the macroscale, we can update the flow and stress fields accordingly. Later in the paper we will provide an example of a lithospheric drip facilitated by lateral compression (e.g. subduction initiation in a passive margin on a plate only undergoing ridge push, as in Fig. 1). Our development here assumes a 2-D domain, although a 3-D development is also feasible and outlined in Appendix  A. However, it is important to reiterate that at the macroscopic scale, the velocity |${\mathbf {v}}$|⁠, strain-rate |$\dot{\underline{\mathbf {e}}}$| and deviatoric stress |${\underline{{\boldsymbol{\tau }}}}$| are treated as if they exist in an incompressible single-phase, albeit effectively anisotropic, medium. Information about separate phases is only retained at the meso- and microscales and, as discussed later, the stresses within separate phases and between mixed and unmixed zones can differ, but these ‘subscale’ stresses are only used to evolve micostructure and infer viscosities at these scales.

2.2.1 The macroscopic stress field as seen by the mesoscale patch

While processes inside the mesoscale patch depend on the macroscale velocity |${\mathbf {v}}$| and its gradients, the details of these kinematic fields are not important yet. However, some generic features of the stress field |${\underline{{\boldsymbol{\tau }}}}$|⁠, specifically its principal values and directions, are important for articulating the physics of the mesoscale patch.

A two-dimensional deviatoric stress field can be expressed as

where |$\hat{{\mathbf {x}}}$| and |$\hat{{\mathbf {z}}}$| are unit vectors, and |$\tau _{xx}$| and |$\tau _{xz}$| are the normal and shear stresses, respectively. The two principal stresses or eigenvalues of this tensor are

where |$\tau ^2$| is the 2nd invariant of the stress tensor. The eigenvectors or principal directions for the stress tensor are defined by the mutually orthogonal unit vectors associated with the compressive (⁠|$-\tau$|⁠) and tensile (⁠|$+\tau$|⁠) eigenvalues, respectively:

where |$-\mathrm{sign}(\tau _{xx})$| is chosen to break the sign degeneracy for eigenvectors, such that |$\hat{{\mathbf {n}}}$| will be positive and point in the compressive stress direction (e.g. if |$\tau _{xx}=-\tau$| and |$\tau _{xz}=0$|⁠). The compressive component of the stress tensor is a dyad composed of the eigenvector |$\hat{{\mathbf {n}}}$| and its associated eigenvalue |$-\tau$|

These features of the stress field are imposed on the mesoscale patch to be described next.

2.2.2 The mesoscale patch

The macroscopic domain is comprised of much smaller units of space that are big enough to include a significant number of parcels of two different phases, and in which each phase is comprised of multiple grains. Since we treat this system in 2-D, we refer to this mesoscale unit of space as a patch.

The mesoscale patch has a size (diameter or edge length) |$\delta$| that is, by definition, much smaller than the length scale of the macroscopic domain, but is also based on the characteristic wavelength |$\lambda$| of the phase heterogeneity. Specifically, |$\delta\ge \lambda$| wherein the patch contains at least one monophase unit, but ideally even more, such that the heterogeneity is roughly periodic across the patch (see Charalambakis 2010). However, while a larger patch captures more phase heterogeneity, it cannot be so large that macroscopic stress and velocity fields vary significantly across it. In principle, and from a practical perspective, a macroscopic grid has grid-spacing |$\delta$| and thus each gridpoint in the macroscopic domain effectively represents one mesoscale patch. Whether one chooses |$\delta=\lambda$| or |$10\lambda$| or |$100\lambda$|⁠, etc., depends on both (a) how well the phase heterogeneity is represented as a simple periodic structure within the patch, and (b) how well the grid resolves macroscopic quantities.

2.2.3 Phase mixing at the mesoscale

The mesoscale patch is comprised of two phases with volume fractions |$\phi _1$| and |$\phi _2$|⁠, or |$\phi =\phi _1$| and |$1-\phi =\phi _2$|⁠. (Phase 1 with volume fraction |$\phi$| is nominally identified as the secondary phase (e.g. pyroxene in the upper mantle) and phase 2 with fraction |$1-\phi$| is the primary phase (e.g. olivine); however either phase fraction can reach 0 or 1 at any given point in space, thus these descriptions are qualitative at best.) The theory of diffusion between phases from grain-mixing was described previously by Bercovici & Mulyukova (2018, 2021) and Bercovici et al. (2023) and states that mixing occurs across phase fronts and along the compressive stress direction:

where |$\mathcal {K}$| is a diffusion coefficient (⁠|$\mathcal {K}\tau$| has units of diffusivity |$\rm m^2\,s^{-1}$|⁠), and can be a function of average local grain size as well as |$\phi$| itself (see most recently Bercovici et al. 2023, Section 2.1 for a complete summary). However, for now we assume that |$\mathcal {K}$| is dominated by its average over the patch domain, and thus is treated as uniform across the patch, depending on its relation to grain size and phase fraction.

2.2.4 The initial phase field

We seek a simple, analytic solution to eq. (5), so that it can be quickly evolved for every mesoscopic patch, and the ease of solution depends partially on the initial condition for |$\phi$|⁠. There is no unique or implicitly ‘correct’ initial field, and the choice is subjective other than seeking the simplest plausible structure that allows a sufficiently fast solution to eq. (5). For now, we assume that since the patch size |$\delta$| is dictated by the characteristic phase heterogeneity wavelength |$\lambda$|⁠, then, within the patch, no other wavelengths are dominant, and all variations in |$\phi$| have wavelength |$\lambda$|⁠. Thus, we express |$\phi$| in terms of 2-D Fourier components each with a wavelength |$\lambda$|⁠, but such that the initial phase field is as isotropic as possible (i.e. there is no unique direction to the initial phase field pattern, or as little direction as possible). Specifically, we write the phase field as a simple cosine Fourier series, but where the wavelength |$\lambda$| is fixed and summation is over the angles of the wave vector |${\mathbf {k}}$|⁠:

where |$\phi ^{\prime }$| is the perturbation against the mean background value |$\phi _0$|⁠, and |$\tilde{{\mathbf {x}}}$| is the position relative to the centre of the patch (see Section 2.2.5 below). Initially, before shearing and diffusion take place, all the |$|\breve{\phi }_j|$| are the same, and the wave-vectors are at evenly spaced angles and obey

where |$k=2\pi /\lambda$|⁠. If |$N=1$|⁠, then |$\phi ^{\prime }$| is a plane wave, which is intrinsically anisotropic at the macroscale and not useful; for |$N=2$|⁠, |$\phi ^{\prime }$| has a square checkerboard structure, which is also not isotropic enough; for |$N=3$|⁠, |$\phi ^{\prime }$| has a hexagonal tile structure, which is reasonably isotropic; for |$N\gt 3$|⁠, |$\phi ^{\prime }$| becomes increasingly isotropic, but this requires retaining more variables for each patch. We find |$N=3$| (or even 4) reaches the right balance of providing a sufficiently isotropic phase field, with nearly circular phase anomalies, while carrying a minimum number of variables (for example, see the phase fields at |$t=0$| in Fig. 3). In general, for this simple wave pattern, requiring that |$\phi$| extends across the full range of |$0\le \phi \le 1$| demands that initially |$\phi _0= N|\breve{\phi }_j|={1}/{_2}$|⁠. Moreover, if eq. (5) is assumed linear with constant coefficients, then |$\phi _0$| remains fixed. The initial phase field eq. (6) is provided as an example, and other choices of the initial |$\phi$| can be assumed, for example a coarse random field as in Bercovici & Mulyukova (2021) and Bercovici et al. (2023), or even from laboratory data or field hand samples. However, expressing |$\phi$| as a simple superposition of a few wave-like structures facilitates the rapidity of the solution, as shown below.

2.2.5 Coordinate transformations for diffusive mixing

To get an analytic solution for eq. (5) we first go into the reference frame of the patch moving at the mean velocity |$\overline{{\mathbf {v}}}$| in the macroscopic reference frame, where the areal average of any quantity Q on a patch is

Say the macroscopic position of any point in a fixed macroscopic frame is |${\mathbf {x}}$|⁠, and the position of the centre of a specific patch is |${\mathbf {X}}$|⁠; in the frame of reference of the patch, these positions are |${\mathbf {x}}-\overline{{\mathbf {v}}}t$| and |${\mathbf {X}}-\overline{{\mathbf {v}}}t$|⁠, respectively. Moreover, the position of a point inside the patch relative to its centre is the difference |$\tilde{{\mathbf {x}}}={\mathbf {x}}-{\mathbf {X}}$|⁠, which is true regardless of frame of reference [note |$\tilde{{\mathbf {x}}}$| was already used above in Section 2.2.4, particularly in (6)]. This simple Galilean transformation to the moving frame is valid within the patch domain of size |$\delta$| for two reasons. First, |$\overline{{\mathbf {v}}}$| is uniform across that domain. Second, we can assume that |$\overline{{\mathbf {v}}}$| is equivalent to the velocity of the centre of the patch, which is the velocity |${\mathbf {v}}$| at that point in macroscopic space that locates the patch; thus, with a creeping flow approximation at the macroscale (see Section 2.5.1), instantaneous acceleration of the patch frame is zero (i.e. |$\mathrm{d}\overline{{\mathbf {v}}}/\mathrm{d}t=0$|⁠). With this first transformation, we are in fact solving the evolution of each patch in a Lagrangian frame.

We then also rotate the reference frame into the principal direction frame of the imposed stress tensor |${\underline{{\boldsymbol{\tau }}}}$|⁠, that is with a unit direction |$\hat{{\mathbf {n}}}$| in the compressive stress direction, and the orthogonal one |$\hat{{\mathbf {m}}}$| in the tensile direction (as described in eq. 3) to give

where |$\tilde{{\mathbf {v}}} = {\mathbf {v}}-\bar{{\mathbf {v}}}$| is the velocity of material in the patch relative to the centre of the (moving) patch, but also transformed into the principal stress coordinate system. The principal stress coordinates are |$x_n$| and |$x_m$| in which the position vector inside the patch relative to the patch’s centre is |$\tilde{{\mathbf {x}}} = \tilde{x}\hat{{\mathbf {x}}}+ \tilde{z}\hat{{\mathbf {z}}}= x_n\hat{{\mathbf {n}}}+ x_m\hat{{\mathbf {m}}}$|⁠. Note that in the transformed coordinate system |${\underline{{\boldsymbol{\tau }}}}_c = -\tau \hat{{\mathbf {n}}}\hat{{\mathbf {n}}}$|⁠, in which |$\tau$| as well as |$\mathcal {K}$| are constant over the patch, and thus the diffusion term on the right-hand side of eq. (5) collapses to that on the right-hand side of eq. (9).

The evaluation of |$\tilde{{\mathbf {v}}}$| involves describing the flow and deformation in the patch, relative to the patch centre, and after transformation to the principal stress frame, which, if done exactly, is cumbersome. Instead, we assume that, since the patch is small, the velocity is well described by the linear expansion

where

are the strain-rate and vorticity tensors, respectively, evaluated from the macroscopic velocity field at the centre of the patch. However, in the coordinate system rotated to the principal stress directions, the strain-rate is defined by pure shear and is |$\dot{\underline{\mathbf {e}}}=-\dot{e}(\hat{{\mathbf {n}}}\hat{{\mathbf {n}}}-\hat{{\mathbf {m}}}\hat{{\mathbf {m}}})$|⁠, while |$\underline{{\boldsymbol{\Omega }}} = \Omega (\hat{{\mathbf {n}}}\hat{{\mathbf {m}}}-\hat{{\mathbf {m}}}\hat{{\mathbf {n}}})$|⁠, where |$\dot{e}^2$| is the second invariant of |$\dot{\underline{\mathbf {e}}}$|⁠, and |$\Omega = -{1}/{_2}\hat{{\mathbf {y}}}\cdot \boldsymbol{\nabla }\times {\mathbf {v}}$| is proportional to the 2-D vorticity; note that |$\underline{{\boldsymbol{\Omega }}}$| has that same form, and |$\Omega$| the same value, for any transformation between 2-D coordinate systems. In this case, the velocity within the 2-D patch from eq. (10) becomes

Also note that the strain-rate tensor is in its principal form when the system is rotated into the principal stress directions provided that the material within the patch is isotropic, albeit heterogeneous; it is only when the heterogeneous structure is averaged over the patch that it has an effective anisotropic fabric, as discussed below in Section 2.2.7.

2.2.6 Analytic solution for the phase field

In addition to the phase-field evolution eq. (9), the Fourier component wave vector |${\mathbf {k}}_j$| also evolves according to conservation of wave fronts given by

(Olson et al. 1984; Whitham 2011) where we drop the subscript j for convenience (for now); with eq. (12), this yields

where |${\mathbf {k}} = k_n\hat{{\mathbf {n}}}+k_m\hat{{\mathbf {m}}}$|⁠. After some tedious math, the solution to eq. (14) yields

where |${\mathbf {K}}$| is the value of |${\mathbf {k}}$| at |$t=t_0$| (e.g. as in eq. 7 for |$t_0=0$|⁠), |$\alpha = \sqrt{\dot{e}^2-\Omega ^2}$| and

This solution, however, is not valid for the special case of |$\dot{e}^2=\Omega ^2$|⁠, as in the example of simple shear; in this situation, the solution to eq. (14) is simply

(Olson et al. 1984). The diffusion eq. (9) for each Fourier mode then reduces to

for which there is an analytic solution

for |$\alpha \ne 0$| and

for |$\alpha = 0$|⁠.

The full phase field within the patch thus evolves according to

where |${\mathbf {k}}_j$| is the same as |${\mathbf {k}}$| above but restoring the subscript j. Since |${\mathbf {k}}_j\cdot \tilde{{\mathbf {x}}}$| is a scalar it does not matter in what coordinate system the vectors in the product are evaluated. Three idealized examples of how |$\phi$| evolves under simple macroscopic stress conditions are shown in Fig. 3. Moreover, note that the above solutions assume that the flow field, which governs |$\dot{e}$| and |$\Omega$|⁠, is constant over the time increment from |$t_0$| to t; in practice the above solutions are applied from one numerical time step of size |$\mathrm{d}t$| to the next, where |$t-t_0=\mathrm{d}t$| is very small (and determined by the CFL condition).

In the case where viscosity variability is only due to one phase being weaker than the other phase (and mixing or weakening associated with mixing is negligible), as in the study by Girard et al. (2016), then the above equations can also be used, but in the limit of |$\mathcal {K}=0$| such that |$\breve{\phi }_j$| remain unchanged, but the wavenumber |${\mathbf {k}}$| still evolves through distortion by the flow field according to eq. (15) or (17).

2.2.7 Fabric tensor

Given an analytic solution for |$\phi$| in each patch, we seek the fabric associated with weak zones, which coincide with bands of mixing between phases wherein grain size is diminished and viscosity reduced (see Section 2.3). The fabric of simple two phase mixtures was defined by Bercovici et al. (2001a) as a tensor based on the orientation of the interface between phases, or phase fronts, given by, for example, the phase-front normal |$\boldsymbol{\nabla }\phi /|\boldsymbol{\nabla }\phi |$|⁠. Specifically, any directional tensorial properties of the two-phase or porous medium, such as anisotropic permeability, are necessarily dictated by the interface orientation (since there is no other direction in the mixture, assuming the phases themselves are materially isotropic), and thus proportional to a fabric tensor. In the simple mixture case, the fabric tensor would go as the dyad |$(\boldsymbol{\nabla }\phi )(\boldsymbol{\nabla }\phi )$|⁠.

The mixing fronts for diffusive grain mixing are, however, normal to the mixing flux vector |$\mathcal {K}{\underline{{\boldsymbol{\tau }}}}_c\cdot \boldsymbol{\nabla }\phi$| (note that mixing fronts only occur where this vector is non-zero) and thus the fabric tensor of the mixing zones is constructed from a dyad composed of this vector:

which is designed so that the trace of |$\underline{\mathbf {F}}^{\prime }$| is also the rate of work done (per unit volume) by diffusive grain mixing (see Bercovici & Mulyukova 2018, 2021; Bercovici et al. 2023). (A somewhat different fabric tensor was defined in Bercovici & Mulyukova (2021) that involved gradients of an isotropic but heterogeneous viscosity field; however, here we neither know a priori, nor want to calculate explicitly, the full isotropic heterogeneous viscosity field in the patch; instead we seek a simpler and faster prediction of the fabric tensor based purely on the analytic solution for |$\phi$|⁠.) The spatial average of the fabric tensor over the patch space yields the fabric for that point in macroscopic space:

and the factor of |${1}/{_2}$| in (23) arises from integration of squared sines (assuming |$\phi ^{\prime }$| is periodic in the patch domain). The eigenvectors of |$\underline{\mathbf {F}}$| give the orientation of the fabric; however, we immediately recognize that these are simply |$\hat{{\mathbf {n}}}$| and |$\hat{{\mathbf {m}}}$|⁠, which are associated with the eigenvalues representing the maximum and minimum mixing flux, |$\mathcal {F}_n=\tfrac{1}{2}\mathcal {K}\tau ^2 \sum _j^{N}k_{nj}^2\breve{\phi }_j^2$| and |$\mathcal {F}_m=0$|⁠, respectively. Thus, no further calculation is necessary to infer that |$\hat{{\mathbf {n}}}$| is the maximum mixing flux direction and is thus the unit normal to the weak zone laminae.

The above fabric tensor is predicated on the assumption that weak zones coincide with mixing zones between phases. In cases wherein the phases have very different viscosities and mixing is negligible (or weakening associated with mixing is less important) then the fabric tensor, as in the simple two-phase mixture (Bercovici et al. 2001a), only depends on the morphology of the interface between phases, and is, again, expressed by the patch-averaged dyad |$\overline{\boldsymbol{\nabla }\phi \boldsymbol{\nabla }\phi } =\tfrac{1}{2} \sum _j^{N}{\mathbf {k}}_j{\mathbf {k}}_j\breve{\phi }_j^2$|⁠. In this case, the eigenvectors and eigenvalues are not as straightforward as in the grain-mixing case above, but are still readily calculated, and one can infer that the weak zones are predominantly aligned normal to the eigenvector associated with the maximum fabric-tensor eigenvalue. In cases where both phase mixing and a strong biphasic viscosity contrast occur, the fabric tensor should be a hybrid or average of these two constructions, weighted according to the respective viscosity contrasts (i.e. the contrast between mixed and unmixed zones versus the contrast between phases).

2.2.8 Volume fraction of mixed zones

The degree of anisotropy in the viscosity field of a mesoscale patch in space depends not only on the orientation of the weak bands, but also on the viscosity variability between weak/mixed and strong/unmixed regions (which we will treat when we step down to the grain-scale space in Section 2.3), and the volume fraction of weak/mixed bands.

The weak/mixed zones only exist where |$\phi (1-\phi )\ne 0$|⁠, and thus we can define the local mixed-zone volume fraction as

which has a maximum value of |$\Phi ^{\prime }=1$| at |$\phi ={1}/{_2}$|⁠, and the exponent |$\lambda$| allows some control on the shape of the function (e.g. |$\lambda\ll 1$| allows a more box-car shape, while |$\lambda\gg 1$| gives a narrow shape centred on |$\phi ={1}/{_2}$|⁠). A mean volume fraction of mixed/weak zones across a patch is, as defined in eq. (8),

For simplicity, given our assumed initial phase field (eq. 6), we assume that |$\phi _0={1}/{_2}$| and the initial extrema of |$\phi ^{\prime }=\pm {1}/{_2}$| such that at |$t=0$|⁠, |$\breve{\phi }_j={1}/{_(2N)}$|⁠. In this case, eq. (25) becomes

The linearization of |$\Phi ^{\prime }$| in eq. (26) is a reasonable approximation inside the mixed zones where |$|\phi ^{\prime }|\ll {1}/{_2}$|⁠, but less so where |$|\phi ^{\prime }|\approx {1}/{_2}$| far from the mixed zones, where |$\phi =0$| or 1. Initially |$\breve{\phi }_j = {1}/{_(2N)}$|⁠, but |$\breve{\phi }_j\rightarrow 0$| as |$t\rightarrow \infty$| (see eq. 19). Thus, the mean mixed zone volume fraction goes from |$\Phi = 1-\lambda/(2N)$| at |$t=0$| to |$\Phi =1$| as |$t\rightarrow \infty$|⁠. We can constrain |$\Lambda$| by assuming that the monophase units are initially completely separated and unmixed such that |$\Phi =0$| at |$t=0$|⁠, in which case we select |$\lambda=2N$|⁠. Of course the treatment of |$\Phi$| is slightly more involved if |$\phi _0\ne {1}/{_2}$| and if we do not linearize |$\Phi ^{\prime }$|⁠.

In cases where phase viscosity contrast is more important than weakening associated with mixing, then the average volume fraction of weak zones over the patch is simply either |$\phi _0$| or |$1-\phi _0$|⁠, depending on which phase is intrinsically weaker.

2.2.9 Summarizing where we stand before we (temporarily) leave mesoscale space

So far we have defined a phase fraction field |$\phi ({\mathbf {x}},t)$| in a small patch of mesoscale space with a simple approximation eq. (6), which can be represented with just a few numbers that we need to evolve with time (i.e. the phase fraction Fourier amplitudes |$\breve{\phi }_j$| and the Fourier wave numbers |${\mathbf {k}}_j$|⁠). With a proper rotation, diffusive grain mixing of phases in the presence of deformation has an analytic solution (e.g. eqs 17, 19 and 21), which allows for the rapid calculation of mixing at each patch (and recall each patch is a point in macro-space). The solution for the phase fraction can also be used to quickly calculate the orientation of the fabric tensor |$\underline{\mathbf {F}}$| from eq. (23), as well as the mixed zone volume fraction |$\Phi$| from eq. (26). These latter two quantities will let us build the anisotropic viscosity field, but we first must infer the viscosity within the mixed and unmixed regions, which also involves determining the grain size evolution in these regions, albeit with rapid solutions. Thus we next venture down another scale into the microscopic grain-space, but we will return later to the mesoscale space to determine the anisotropic viscosity field.

2.3 From the mesoscale patch to microscale grain space

There are three types of regions in the mesoscale patch. There are monophase zones, or tectonites, of each type of phase, that is with |$\phi =0$| or |$\phi =1$|⁠, and there are mixed zones between the tectonites with |$0\lt \phi \lt 1$|⁠. We assume that, within a given patch, each type of region undergoes the identical evolution of grains and hence viscosity (and in principle the patch is all at one temperature).

2.3.1 Mixed zones undergo two-phase grain damage

Grain size evolution in the mixed zones, where the two phases coexist, is assumed to follow two phase grain damage theory (which describes the competition between healing via surface-tension driven coarsening, and deformationally driven damage and grain comminution), enhanced by Zener pinning (Bercovici & Ricard 2012, 2013, 2014, 2016). This theory was most recently adapted to include diffusive mixing by Bercovici & Mulyukova (2018, 2021) and Bercovici et al. (2023) and appears as

where |$R_i$| and r are the grain sizes of phase i and the interface roughness (radius of curvature); |$G_i$| and p are the rate and the exponent of grain boundary coarsening, respectively, and |$G_I$| and q are the analogous values for interface coarsening; |$\mathcal {Z}_i$| is the Zener pinning factor in which |$\mathfrak {c}= 0.87$| [associated with the assumed grain size distribution; see Bercovici & Ricard 2012, eq (8)]; |$\gamma _I$| and |$\gamma _i$| are the interface and grain boundary surface energies, respectively; and |$\mathfrak {n}=3\phi (1-\phi )$| represents the dependence of interface area density on volume fraction (based on a spherical pore model in which the area density is approximately |$\mathfrak {n}/r$|⁠). The partitioning factors |$f_I$| and |$f_G$| determine the fraction of deformational work going into new interface and grain boundary area, respectively. The work to produce more grain boundary area occurs via dynamic recrystallization and is restricted to dislocation creep work |$2A_i\tau _i^{n+1}$|⁠, where |$A_i$| and n are the dislocation-creep compliance and power-law index for phase i, and |$\tau _i$| is the local invariant stress (square root of the 2nd stress invariant) of phase i in the mixed zones. The work to produce more interface area is extracted from the total work rate, that is the bulk deformational work |$2\tau _{\upsilon }\dot{e}$| (where |$\tau _{\upsilon }=\sum _i\phi _i\tau _i$| is the volume averaged local invariant stress) and that associated with grain-scale mixing, since deformation by any creep mechanism will cause interface distortion by stretching, rending or mixing. The mixing work is assigned to be equivalent to |$2\mathcal {F}_n$| (see eq. 23), where the factor of 2 accounts for |$\mathcal {F}_n$| being the patch-average, not the peak mixing work expected in the mixed zones. Although the interface surface energy is proportional to |$\gamma _I\mathfrak {n}$|⁠, we include a factor in the interface damage term |$\mathfrak {n}^{\ell +1}$| adjacent to |$f_I$| (in which |$\ell \gt 0$|⁠) since we expect no interface damage when |$\phi =0$| or 1 hence when |$\mathfrak {n}=0$|⁠. As in our previous studies (most recently Bercovici & Mulyukova 2021; Bercovici et al. 2023), the factor |$\zeta$| equals 4.95, but is specific to the assumed lognormal grain size distribution.

These equations are posed in the frame of reference of the patch and thus the advection velocity is approximately 0 (although the local velocity is that of uniform pure shear and rotation, the average velocity is zero). Note also that the local stresses |$\tau _i$| and |$\tau _{\upsilon }=\sum _i\phi _i\tau _i$| used in the above equations are associated with the grain size in the mixed zones through a composite constitutive law for a given strain-rate |$\dot{e}$|⁠, as discussed below in Section 2.3.3. The phase fraction in the mixed zone is |$\phi \approx {1}/{_2}$| (not necessarily |$\phi _0$|⁠, since locally the mixed zones only see boundary conditions of |$\phi =0$| and |$\phi =1$|⁠, not the bulk composition, although in this paper we have assumed the bulk composition has |$\phi _0={1}/{_2}$| anyway). For the mixed zone, the output of the evolution model are the mixed-zone grain sizes |$R_{\sf {m}i}= R_i$| obeying eq. (27).

2.3.2 The monophase zones undergo single-phase grain damage

Grain size evolution in the mono-mineralic zones with |$\phi =0$| or 1 is governed by single-phase damage (Ricard & Bercovici 2009; Rozel et al. 2011), which is similar to the palaeowattmeter relations of Austin & Evans (2007). Since grain-damage is entirely dependent on dynamic recrystallization, grain size reduction is restricted to the dislocation creep regime, and grains tend to be larger and the mono-phase zones stronger. In this case, grain size evolution obeys

As with the mixed zone, the local stress |$\tau _i$| used in the above equation is associated with the grain size in the unmixed zones through a composite constitutive law for a given strain-rate |$\dot{e}$| (see Section 2.3.3). For the unmixed zone, the unmixed-zone grain sizes are identified as |$R_{\sf {u}i}= R_i$| obeying eq. (28).

2.3.3 Viscosity inside and outside the mixed zones

To be consistent with the patch evolution model that permits an analytical solution (Section 2.2.3), we assume the invariant strain-rate |$\dot{e}$| is, at the small mesoscale, uniform across the patch, and thus the submesoscale mixed and unmixed zones experience the same deformation. This assumption imposes strain compatibility within a mesoscale patch, such that the weak and strong zones cannot separate and open voids and gaps. But the assumption also means the mechanical response of submesoscale weak zones is expressed in terms of reduced stress (rather than being deformed faster than the adjacent strong zones). However, while homogeneous strain in the patch disallows strain localization on weak mixed zones at the subpatch scale, it permits stress and thus viscosity reduction in these zones, which thus still resolves weakening fabric and the aggregate loss of strength of the patch. Moreover, deformation is not uniform at the macroscale, that is it can still vary across multiple patches (just not within a patch), which is consistent with the continuum approximation. Thus if a patch or group of adjacent patches weakens in bulk (i.e. because of the summative effect of the weak-zones within the patch) then the net effective viscosity (which as shown can be anisotropic) is reduced. If the system is driven by a buoyancy or stress load, the weakened zone’s response would be to increase the patch’s bulk strain-rate, which, at the macroscale might appear as a shear-localized zone resolved across multiple patches. In the same sense, the evolution of microstructure within a patch needs to be continuous between patches, in which case the strain-rate imposed on adjacent patches needs to be continuous as well. In short, the patches cannot be so big that they cannot resolve gradients in velocity, which again is simply specifying that a patch is the equivalent to a canonical continuum-mechanical control volume.

2.3.4 Viscosity inside the mixed zones

The mean effective viscosity within a mixed zone depends on the stress |$\tau _i$| within that weak zone for phase i, which is the solution to the composite constitutive law for diffusion and dislocation creep

where |$B_i$| and m are the diffusion creep compliance and grain size exponent, respectively. Eq. (29) can be quickly determined analytically (for |$n=1,2$| and 3) or numerically with a rapid approximation scheme (see Bercovici et al. 2023, Appendix A1.1). Since the phases are well dispersed with each other inside the mixed zone, the mean viscosity in that zone is

where again we have assumed the phase fraction of the mixed zone is on average |$\phi ={1}/{_2}$|⁠.

2.3.5 Viscosity in the unmixed zones

The viscosity of an unmixed zone likewise depends on the stress in that zone, which is the solution to

such that viscosity of the unmixed phase is

The mean viscosity of the unmixed zones averaged over the bulk of the patch is

which for |$\phi _0={1}/{_2}$| yields the same relation as eq. (30) except for the stress associated with the larger grain sizes |$R_{\sf {u}i}$|⁠.

However, for the problem in which one phase is weaker than the other, and mixing (or weakening associated with mixing) is negligible, then we would preserve both phase viscosities |$\mu _i$| for later determination of mechanical anisotropy.

Finally we reiterate that the stresses calculated for these viscosities are also used in the mixed and unmixed zone grain-evolution laws eqs (27) and (28).

2.4 Back to mesoscale space to infer an anisotropic viscosity field

Given the estimates of viscosities in the weak and strong regions determined at the grain scale, we return to the mesoscale patch to infer an anisotropic viscosity.

We assume that the weak bands resulting from both diffusive mixing (with stress dependent anisotropic diffusivity) and shearing leads to an overall bilaminated effective structure (see Fig. 3), in which the effective viscosity for compressive or tensile deformation normal to the laminae, and that associated with shear deformation parallel to the laminae are (Treagus 2003)

respectively. For the case in which mixing is neglected and the phases have different viscosities, we would replace |$\Phi$| with |$\phi _0$|⁠, |$\mu _{\sf {m}}$| with |$\mu _1$|⁠, and |$\mu _{\sf {u}}$| with |$\mu _2$| in (34). Other non-laminate structures, such as ellipsoidal phase inclusions, can also be considered (Treagus 2002, 2003; Fletcher 2004; Montési 2007; Dabrowski et al. 2012), which are considerably more complex and require evolution of other parameters, namely the ellipticity (see discussion by Dabrowski et al. 2012). However, elliptical or ellipsoidal axes most likely follow the strain principal directions, and thus reach high ellipticities and nearly laminate structures by strains |$\gt 1$|⁠, and even faster if the weak-zones coincide with mixing regions, which begin semi-laminated (see Fig. 3). Nevertheless, treating anisotropy with phase structures that are non-laminated (elliptical/ellipsoidal) would be a worthy expansion of this step in our ‘upscaling’ framework

The effective laminae are assumed to be oriented according to the fabric tensor |$\underline{\mathbf {F}}$|⁠, but this simply means that the laminae are normal to the principal compressive stress direction |$\hat{{\mathbf {n}}}$| and parallel to the conjugate vector |$\hat{{\mathbf {m}}}$|⁠. (In principle, the magnitude of the fabric tensor |$\mathcal {F}_n$| also indicates the degree of lamination; however it is operationally redundant with mixing volume fraction |$\Phi$| because it also indicates the presence of mixing zones. For example, |$\mathcal {F}_n = 0$| implies there is no diffusive mixing between phases, which is also indicated by |$\Phi =0$|⁠.) Thus the normal viscosity |$\mu _{\rm{\small {N}}}$| is the effective viscosity for deformation under normal stress in either the |$\pm \hat{{\mathbf {n}}}$| or |$\pm \hat{{\mathbf {m}}}$| direction (i.e. compression or tension is limited by the response of the stronger medium, depending on its volume fraction), and |$\mu _{\rm{\small {S}}}$| is the effective viscosity for deformation under shear parallel to the laminae (i.e. tractions acting on surfaces pointing in the |$\pm \hat{{\mathbf {n}}}$| direction, but dragging in the |$\pm \hat{{\mathbf {m}}}$| direction).

The normal and shear stress components are generally

where |${\underline{{\boldsymbol{\tau }}}}$| is the macroscopic stress tensor defined in eq. (1). The unit normals |$\hat{{\mathbf {n}}}=a\hat{{\mathbf {x}}}+b\hat{{\mathbf {z}}}$| and |$\hat{{\mathbf {m}}}=-b\hat{{\mathbf {x}}}+a\hat{{\mathbf {z}}}$|⁠, as defined in eq. (3), lead to

The normal and shear strain-rate components of strain-rate |$\dot{\underline{\mathbf {e}}}$| obey a similar transformation.

The bilaminate structure has a mechanical response that obeys a transverse isotropic rheology, dictated by the normal and shear viscosities |$\mu _{\rm{\small {N}}}$| and |$\mu _{\rm{\small {S}}}$| (see, e.g. Honda 1986):

where |$\dot{e}_{nn}$| and |$\dot{e}_{nm}$| are the normal and shear strain rates relative to the laminae orientations. Restating eq. (37) in terms of Cartesian components

and converting to the original Cartesian coordinates, we arrive at an anisotropic rheological law

where the associated anisotropic viscosity due to the bilaminated mixed layer structures is

where |$\Delta \mu = \mu _{\rm{\small {N}}}-\mu _{\rm{\small {S}}}$|⁠. The relation for |$\underline{\mathbf {M}}$| can be written more compactly, while also separating its isotropic and anisotropic components:

where |$\bar{\mu }=(\mu _{\rm{\small {N}}}+\mu _{\rm{\small {S}}})/2$| and

It is worth noting a few important limits implied by eq. (41). At time |$t=0$|⁠, before mixing begins, |$\Phi \approx 0$| in which case |$\mu _{\rm{\small {N}}}=\mu _{\rm{\small {S}}}=\bar{\mu }=\mu _{\sf {u}}$| (see (34)), and thus |$\Delta \mu =0$|⁠, such that the viscosity is isotropic and at the large value of the strong, unmixed grains. As time progresses, and |$\Phi \rightarrow {1}/{_2}$|⁠, then |$\bar{\mu }$| decreases and |$\Delta \mu$| increases, and the anisotropy is maximized (i.e. the viscosity ratio |$\mu _{\rm{\small {N}}}/\mu _{\rm{\small {S}}}$| is maximized at |$\Phi ={1}/{_2}$|⁠). As |$t\rightarrow \infty$|⁠, |$\Phi \rightarrow 1$|⁠, whereupon |$\mu _{\rm{\small {N}}}\approx \mu _{\rm{\small {S}}}\approx \bar{\mu }\rightarrow \mu _{\sf {m}}$|⁠, and |$\Delta \mu \rightarrow 0$|⁠; the patch becomes isotropic again, but its viscosity is that of the weak, well-mixed, small-grain sized medium, |$\mu _{\sf {m}}$|⁠, as expected. Thus, in short, before mixing starts, the medium initiates isotropically with a high viscosity |$\mu _{\sf {u}}$|⁠, and once mixing is complete it ends isotropically with a low viscosity |$\mu _{\sf {m}}$|⁠; and at intermediate times it is anisotropic.

For cases in which the phase viscosities are comparable, we can further simplify the expression for |$\underline{\mathbf {M}}$| for the sake of mathematical compactness and pose the viscosity field in terms of a reference viscosity. Specifically, we assume there is a reference viscosity |$\mu _0$| that exists before mixing begins, and is associated with a given grain size |$R_0$| and stress |$\tau _0$| such that

where the far right-side approximation is made assuming the two phase rheologies are similar. The reference grain size is assumed to be the steady-state solution for mono-phase grain damage in each phase (equivalent to the piezometric or palaeowattmetric grain size) from eq. (28)

where the right side approximation is made assuming the properties of the phases do not differ significantly. (In the case when the phase rheologies and grain-growth coefficients differ significantly, we can also express eqs (43) and (44) in terms of phase-averaged properties, which is slightly more cumbersome.) The isotropic viscosity |$\bar{\mu }$| can be expressed relative to this reference viscosity as

and we similarly define |$\Delta \nu =\Delta \mu /\mu _0$|⁠. Further noting that |$\Delta \nu$| and |$\underline{\mathbf {J}}$| never appear separately, we can rewrite eq. (41) as

is the Gary tensor, which uniquely defines the anisotropy of the viscosity field. Note that |$\nu$| and |$\Delta \nu$| represent two important viscosity differences: |$\Delta \nu$| represents the difference in viscosity between laminar layers that lead to viscous anisotropy; |$\nu$| represents the change in isotropic viscosity of the basic state due to mixing and damage at the meso- and microscales.

Finally, we note the implied constitutive law |${\underline{{\boldsymbol{\tau }}}}= 2\underline{\mathbf {M}}\cdot \dot{\underline{\mathbf {e}}}$| is not written as 2nd rank tensor relation, but as an effective 1st rank tensor relation using only the independent components of the stress and strain-rate tensors as pseudo-vectors, as in eqs (37)–(39). Reposing |$\underline{\mathbf {M}}$| as a 4th rank tensor to apply to the full 2nd rank stress and strain-rate tensors is not difficult, just cumbersome and unnecessary. Although so far we have posed the model in 2-D, it is similarly extendable to 3-D (see Appendix  A).

2.4.1 Returning (some) quantities to the Eulerian frame

Most of the quantities such as, for example, the fabric tensor |$\underline{\mathbf {F}}$| eq. (23), mixing volume fraction |$\Phi$| eq. (26) and viscosity |$\underline{\mathbf {M}}$| eq. (46) are instantaneously determined based on other time dependent quantities such as phase fraction |$\phi$| eq. (21), grain sizes |$R_i$| (in the mixed and unmixed zones, namely |$R_{\sf {m}i}$| and |$R_{\sf {u}i}$|⁠), and interface roughness r eqs (27) and (28). The updated values of |$\phi$|⁠, |$R_i$| and r at time t at each patch, or equivalent point in macrospace, are retained and used to initiate the next time step in their evolution. However, they also need to be transformed back to the original macroscopic Eulerian frame from the patch-centred Lagrangian frame, so that they can be updated at the next time step with the new velocity and stress fields. The wave vector |${\mathbf {k}}$| in the relation for |$\phi$| eq. (21) is also transformed back from the principal stress coordinate system to the original macroscopic grid system, so that it can be transformed again to new principal stress directions at a later time. After determining how a scalar, for example |$R_i$|⁠, has evolved in the Lagrangian moving frame of reference from t to |$t+\mathrm{d}t$|⁠, its position is effectively shifted to |${\mathbf {x}}+\overline{{\mathbf {v}}}\mathrm{d}t$|⁠, and similarly for other time-varying quantities. A numerically efficient scheme to execute this step is the method of backward characteristics (e.g. Baptista 1987; Malevsky & Yuen 1991; De Smet et al. 2000). In this approach, the updated Eulerian grid value of, for example |$R_i({\mathbf {x}}, t+\mathrm{d}t)$| is assigned the pre-advected value |$R_i({\mathbf {x}}-\overline{{\mathbf {v}}}\mathrm{d}t, t+\mathrm{d}t)$|⁠, which is determined from the interpolation of |$R_i$| from neighbouring gridpoints. The interpolation (e.g. with a bilinear interpolation) from surrounding gridpoints assumes that |$R_i$|⁠, as well as other time dependent quantities such as |$\phi$| and r, are continuous macroscopic functions at all times (which is justifiable since their evolution is fundamentally being forced by continuous macroscopic stress and velocity fields, provided that their initial conditions are continuous as well). This assumption is also consistent with our prescription that the area of a mesoscopic patch is equivalent to a rectangle defined by four macroscopic gridpoints, although the patch is actually centred on a gridpoint (Section 2.2.2). Other methods that use both Lagrangian and Eulerian formulations can be well suited to model geodynamic flows across multiple scales, such as particle tracking (see van Keken et al. (1997) and an example applied to grain-damage theory by Bellas et al. (2018)), and specifically Lagrangian methods developed for modelling large strains that enable accurate tracking of deformation history (deformable particle-in-cell method, e.g. Moresi et al. 2003; Samuel 2018).

2.5 Back up to macroscopic space to determine stress and flow

Given a relation for the anisotropic viscosity field at every mesoscale patch, and thus equivalently at each point in macroscopic space, we can calculate the macroscopic stress and flow fields.

2.5.1 Governing flow equations

The general flow equations, assuming an incompressible, highly viscous medium, are given by mass and momentum conservation, or equivalently continuity and creeping-flow force balance:

where P, |$\rho$| and g are pressure, density and gravity, respectively, and

2.5.2 Reduction to a single equation for stream-function

For our 2-D model in the x-z plane, we can take |$\hat{{\mathbf {y}}}\cdot \boldsymbol{\nabla }\times$| of eq. (48) to obtain

where the differential operator

differs slightly (by a factor of |${1}/{_2}$|⁠) from its original definition by Bercovici (1993). We have also opted to adhere to strict incompressibility (i.e. not a Boussinesq approximation) by keeping density constant; thermal or chemical buoyancy is easily added by including terms of order |$\partial \rho /\partial x$|⁠, but then a relation for evolving the density (by heat or mass transport) would need to be added.

Using eq. (46) for |$\underline{\mathbf {M}}$| in eq. (49), the stress components are

where |$\Gamma _{ij}=\Delta \nu J_{ij}$| and the |$J_{ij}$| are given by eq. (42). Further satisfying mass conservation eq. (47) by defining |${\mathbf {v}} = \boldsymbol{\nabla }\times (\psi \hat{{\mathbf {y}}})$| in terms of a stream-function |$\psi$|⁠, we can write

and similarly the one component of the vorticity tensor (see eq. 11) is

Using eqs (52) and (53), eq. (50) becomes

which can be reorganized to

Thus given the anisotropic viscosity determined from the meso- and microscales, we infer the macroscopic flow from eq. (56) in terms of a stream function |$\psi$|⁠, with which we infer the macroscopic velocity |${\mathbf {v}}$|⁠, strain-rate |$\dot{\underline{\mathbf {e}}}$|⁠, stress |${\underline{{\boldsymbol{\tau }}}}$| and vorticity tensor |$\underline{{\boldsymbol{\Omega }}}$|⁠; these macroscopic fields are used to update the evolution of |$\phi$| and |$R_i$| at the meso- and microscales, after which the viscosity field is updated, and so on, in a cycle of going down and back up in scales.

Given our upscaling framework to go from macrostructure to microstructure and back again, we provide a working example in the next section.

3 A WORKING EXAMPLE, THE DRIP INSTABILITY WITH LATERAL COMPRESSION: A SIMPLE MODEL OF PASSIVE MARGIN COLLAPSE AND SUBDUCTION INITIATION

3.1 Model set-up

We consider a simple model of ocean lithosphere near a passive margin and the conditions for it to undergo gravitational instabilities leading to subduction (Fig. 1). The model lithosphere is a layer of mean thickness H (wherein |$z\in [-H, 0]$|⁠) and length |$L\gt H$| (wherein |$x\in [-L/2, L/2]$|⁠) that is bounded above by an inviscid atmosphere and below by a relatively inviscid asthenosphere (Fig. 4). The lithosphere has density |$\rho$| that is |$\Delta \rho$| greater than that of the asthenosphere below. The layer is under an imposed, fixed lateral compressive stress |$\tau _0$| representative of ridge push stress (see Dahlen 1981; Mulyukova & Bercovici 2018b; Bercovici & Mulyukova 2021) or intracontinental convergence (e.g. Conrad & Molnar 1997; Molnar & Houseman 2004), and can eventually develop a gravitational Rayleigh–Taylor drip instability on the bottom boundary of thickness h, which represents a downwelling drip, or subduction initiation (of sorts).

The goal of the model is to provide a simple example of the upscaling scheme, but also to compare the growth rate of the drip instability between the simple isotropic case, and the case in which anisotropic weakening has evolved through processes at the meso- and microscales. To keep the demonstration simple we follow a modal analysis in which we assume (or impose) that at the macroscale, the flow field and bottom boundary deflection have a single sinusoidal mode only, for example the deflection |$h = \breve{h}(t)\cos (\boldsymbol{k}x)$|⁠, where for example, |$\boldsymbol{k}= 2\pi /L$| (i.e. one sinusoidal mode across the length L). (Note that |$\boldsymbol{k}$| here is not the same as k in the description of the phase field |$\phi$|⁠.) To preserve the modal symmetry, all the macroscale equations are linearized, around a background state that is isotropic and of uniform grain size |$R_0$|⁠, subject to lateral compressive stress |$\tau _0$|⁠, and thus with an associated isotropic viscosity |$\mu _0$|⁠. Perturbations to this state are caused by the development of weak zones and viscous anisotropy from meso- and microscale processes, and also by the deflection of the bottom boundary. Thus the relative viscosity perturbations that scale as |$\nu$| and |$\Delta \nu$| in eq. (46) (recalling that |$\underline{{\boldsymbol{\Gamma }}}=\Delta \nu \underline{\mathbf {J}}$|⁠) are assumed |$\ll 1$|⁠, and the bottom boundary deflection is |$h\ll H$|⁠.

The protocol for this model follows the general scheme articulated in Section 2, the only extra wrinkle that warrants elaboration is the treatment of the macroscale model, given the specifics of the geometry, forcing, boundary conditions and linearization. The mesoscale and microscale (grain-scale) treatments that eventually lead to an anisotropic viscosity require no extra elaboration, other than non-dimensionalization (see Section 3.2). A table of (dimensional) model properties and parameters is provided in Table 1.

3.1.1 Macroscale stress and viscosity

The background driving force in the model, even before a drip starts, is from an imposed lateral compressive stress field of magnitude |$\tau _0$| (Fig. 4). Both the drip and development of viscosity variations due to grain-damage and mixing induce a first order stress perturbation |${\underline{{\boldsymbol{\tau }}}}^{\prime }$|⁠. The full stress tensor is thus

and the associated strain-rate tensor is

where |$\dot{e}_0$| is the zeroth order strain-rate (from lateral compression) and |$\dot{\underline{\mathbf {e}}}^{\prime }$| is the first order strain-rate perturbation. In this case, eq. (52) leads to |$\dot{e}_0 = \tau _0/(2\mu _0)$| for the zeroth order constitutive relationship, and the first order perturbation constitutive relation is (since |$\nu$| and |$\Gamma _{ij}$| are small perturbation quantities)

where again |$\mu _0$| is given by eq. (43) with eq. (44).

3.1.2 Flow at the macroscale

The velocity is comprised of zeroth order flow associated with pure shear at a rate |$\dot{e}_0$|⁠, and first order perturbation according to

such that

Similarly, the independent strain components from eq. (53) or eq. (58) become

The background horizontally compressive ‘ridge push’ flow is irrotational and thus the vorticity tensor |$\underline{{\boldsymbol{\Omega }}}$| is only affected by the flow perturbation; that is

Since the zeroth order stress field associated with |$\tau _0$| is uniform, we can substitute eq. (59) directly into eq. (50), and given eq. (62)—recalling that |$\mu _0$| is uniform—arrive at

which is, of course, just the linearization of eq. (56). For our modal analysis, we assume |$x=0$| is a symmetry plane and thus |$\psi ^{\prime }=\breve{\psi }(z)\sin (\boldsymbol{k}x)$|⁠, in which case eq. (64) becomes an ordinary differential equation

where the only modes of |$\nu$|⁠, |$\Gamma _{xx}$| and |$\Gamma _{xz}$| that contribute as forcing functions are

and the coefficients |$\breve{\nu }, \breve{\Gamma }_{xx}$| and |$\breve{\Gamma }_{xz}$| can be found with basic sine or cosine transforms (e.g. |$\breve{\nu }= \tfrac{2}{L}\int _{-L/2}^{L/2}\nu \cos (\boldsymbol{k}x)\mathrm{d}x$|⁠). In the end, eq. (65) is a 4th order equation for |$\breve{\psi }$| subject to four boundary conditions at the top and bottom of the layer. These equations can be generalized to include multiple modes (summing over a range of |$\boldsymbol{k}$|⁠), but the single mode suffices for the sake of illustration.

3.1.3 Macroscale boundary conditions

The model lithospheric layer is assumed to be much more viscous than the underlying asthenosphere, and certainly more so than the ocean or atmosphere above; thus the layer is bounded by relatively inviscid media such that its top and bottom boundaries can be treated as free-slip (shear-stress free). Moreover, we assume that the density contrast between the layer and the overlying medium is very large and that deflections of that boundary are negligible; thus the top is also treated as impermeable (zero vertical velocity). The bottom boundary however supports deflections of thickness |$h(x,t)$|⁠, which influence the stress and velocity boundary conditions. The boundary conditions are fully developed and explained in Appendix  B, and summarized below.

At the top boundary, where |$z=0$|⁠, the impermeable condition requires that the vertical velocity is zero, and the free-slip condition requires zero shear stress, which leads to

At the bottom boundary, where |$z=-H-h$|⁠, conditions on vertical velocity (which is related to the rate of the boundary’s deflection), and continuity of shear and normal stresses, are projected onto the mean boundary |$z=-H$| (see Appendix B2), and lead to, respectively,

where again |$(h,\nu ,\Gamma _{xx}) = (\breve{h}, \breve{\nu },\breve{\Gamma }_{xx})\cos (\boldsymbol{k}x)$|⁠, while |$(\psi , \Gamma _{xz}) = (\breve{\psi }, \breve{\Gamma }_{xz})\sin (\boldsymbol{k}x)$|⁠.

3.1.4 A note on principal stress directions

As noted many times already, the macroscopic flow and stress fields are used to evolve the mesoscale heterogeneity, and require the principal stress directions |$\hat{{\mathbf {n}}}$| and |$\hat{{\mathbf {m}}}$| as given in eq. (3), and the components of these unit vectors depend on quantities such as |$\sqrt{\tau \pm \tau _{xx}}$| where |$\tau =\sqrt{\tau _{xx}^2+\tau _{xz}^2}$|⁠. The linearized stress eq. (57) implies that

The full (not linearized) expression for |$\tau$| is

Although we linearize the stress components above, the linearized version of eq. (70), given by |$\tau \approx \tau _0 - \tau ^{\prime }_{xx}$|⁠, fails to resolve the principal stress directions needed in the mesoscale patch. Specifically, the linearized |$\tau$| is equal to |$-\tau _{xx}$|⁠, and thus one of the components of the unit vectors |$\hat{{\mathbf {n}}}$| and |$\hat{{\mathbf {m}}}$| are always zero [i.e. |$b\propto \sqrt{\tau +\tau _{xx}}$| in eq. (3)], and thus these vectors cannot resolve any differences from the Cartesian unit vectors |$\hat{{\mathbf {x}}}$| and |$\hat{{\mathbf {z}}}$|⁠. Moreover, the linearized |$\tau$| cannot satisfy the requirement that |$\tau _{xz}^2 = \tau ^2-\tau _{xx}^2$| if |$\tau _{xz}\ne 0$|⁠. If we keep the full expression for |$\tau$| and take a Taylor expansion of it to order |$\tau ^{\prime 2}/\tau _0^2$|⁠, then |$\tau +\tau _{xx}= {1}/{_2}{\tau ^{\prime 2}}/\tau _0 \ne 0$|⁠; moreover |$\tau ^2-\tau _{xx}^2 = \tau ^{\prime 2} - \tau _{xx}^{\prime 2} = \tau _{xz}^{\prime 2} = \tau _{xz}^2$| as desired. The bottom line is that when inferring the principal stress directions for the mesoscale patch, we use the full non-linear expression eq. (70) for the invariant stress |$\tau$|⁠.

3.2 Dimensionless governing equations

Non-dimensionalization of the system’s governing equations involves rescaling dependent and independent variables by the intrinsic time, length and mass—or, more practically, stress-scales of the system. The macro-, meso- and microscopic levels of our system all follow the same time and stress scale, but each have different intrinsic length scales.

For the sake of simplicity, we assume the creep compliances of the two phases are the same, such that |$A_i=A$| and |$B_i=B$|⁠. (If we consider cases wherein the phases have very different viscosities, this assumption would need to be altered, and we would instead scale with volume averages such as |$(A,B)=\sum _i\phi _i(A_i,B_i)$| (e.g. Bercovici & Ricard 2013, 2014).) For all scales, we non-dimensionalize stress |${\underline{{\boldsymbol{\tau }}}}$| by the basic driving stress |$\tau _0$|⁠, strain-rate |$\dot{\underline{\mathbf {e}}}$| and the vorticity tensor |$\underline{{\boldsymbol{\Omega }}}$| by |$A\tau _0^n$|⁠, and time by |$\mathsf {t_s}= (A\tau _0^n)^{-1}$|⁠. We herein consider non-dimensionalization of the meso-, micro- and macroscale levels in the same order that the upscaling framework has been described. A table of dimensionless variables and properties are given in Table 2.

3.2.1 Dimensionless mesoscale equations

The mesoscale patch is of characteristic size |$\delta$|⁠. As discussed in Section 2.2.2, |$\delta$| is determined by the intrinsic length scale of phase heterogeneity, which is the heterogeneity characteristic wavelength |$\lambda$|⁠, such that the wave number of the initial heterogeneity is |$k=2\pi /\lambda$|⁠. In short, |$\delta\ge \lambda$| for the patch to capture enough heterogeneity for it to be approximately periodic across the patch domain. Thus we can define |$\delta=\mathsf {n}\lambda$| where |$\mathsf {n}\ge 1$| (and typically an integer), in which case |$k=2\pi \mathsf {n}/\delta$|⁠. Non-dimensionalizing mesoscale distance by |$\delta$|⁠, the dimensionless wavenumber |${\mathbf {k}}$| still evolves according to eq. (15) or eq. (17), but where |$\alpha t$| or |$\dot{e}t$| is already dimensionless (regardless of how |$\alpha$|⁠, |$\dot{e}$| or t are scaled) and |${\mathbf {K}}$| is the initial value of |${\mathbf {k}}$| at |$t=t_0$| [as in eq. (7) when |$t_0=0$|], but is dimensionless and with magnitude |$\mathsf {k}=2\pi \mathsf {n}$|⁠. The phase variation amplitudes still follow eq. (19) or (20) but where |$\tau$| is non-dimensionalized by |$\tau _0$|⁠, |$\dot{e}$| and |$\Omega$| by |$A\tau _0^n$|⁠, and thus the dimensionless phase mixing diffusion coefficient is

The dimensionless fabric tensor following eq. (23) is

the trace of which |$\mathcal {F}_n$| is the diffusive mixing work rate, which has been non-dimensionalized by |$A\tau _0^{n+1}$| the same as all other deformational work terms (see Section 3.2.2). The mixing volume fraction |$\Phi$| from eq. (26) is dimensionless already and remains unchanged. Note that none of the quantities that eventually inform the anisotropic structure of the patch retain any information about the patch size |$\delta$| other than how it is absorbed into the dimensionless mixing diffusion coefficient |$\chi$|⁠.

3.2.2 Dimensionless microscale equations

For the microscopic system, the relevant intrinsic length scale is the grain size at the field boundary between diffusion and dislocation creep at the driving stress |$\tau _0$|⁠. Thus grain sizes (including |$R_0$|⁠) and interface roughness are non-dimensionalized by |$\mathsf {R_s}= (B/(A\tau _0^{n-1}))^{1/m}$|⁠. Using the time, and stress scales as well, the dimensionless version of two-phase grain damage eq. (27) in the mixed zones become

where |$\mathcal {C}_i$| and |$\mathcal {C}_I$| are the coarsening numbers for grains and interface, respectively, and |$\mathcal {D}_i$| and |$\mathcal {D}_I$| are the damage numbers for the grains and interface, respectively, defined as

where recall |$\mathsf {t_s}= (A\tau _0^n)^{-1}$|⁠. Moreover, recall that in the mixed zones we assume |$\phi ={1}/{_2}$| and hence |$\mathfrak {n}={3}/{_4}$|⁠. The solution to the two-phase damage equations yields the mixed-zone grain sizes |$R_{\sf {m}i}$|⁠.

For the mono-phase unmixed zones, the dimensionless version of the single-phase grain-damage relation eq. (28) becomes

where the dimensionless numbers remain as defined in eq. (74). The solution to the one-phase damage equations yields the unmixed grain sizes |$R_{\sf {u}i}$| for each phase.

Moreover, the non-dimensional version of the reference grain size |$R_0$| from eq. (44) arises from the steady state solution for eq. (75) for a dimensionless stress of |$\tau _i=1$|⁠, and leads to

if the coarsening and damage numbers are similar for each phase.

The grain sizes in the mixed and unmixed zones are used to infer the viscosity in each phase |$\mu _i=\tau _i/(2\dot{e})$| within those zones by first finding the stresses in each phase |$\tau _i$| that are the solution to the dimensionless version of the constitutive law eq. (29)

where |$R_i$| is either |$R_{\sf {m}i}$| or |$R_{\sf {u}i}$|⁠. The resulting phase viscosities are then used to infer the anisotropic viscosity structure as already articulated in Sections 2.3.3 and 2.4, which is then used to calculate the macroscopic flow and stress fields.

3.2.3 Dimensionless macroscale equations

For the macroscopic level, distance is non-dimensionalized by the layer depth H, and thus velocity by |$HA\tau _0^n$| and stream function by |$H^2A\tau _0^n$|⁠. The domain is now given by |$z\in [-1, 0]$| and |$x\in [-L/(2H), L/(2H)]$|⁠. The dimensionless equations for stress become |${\underline{{\boldsymbol{\tau }}}}_0 = -(\hat{{\mathbf {x}}}\hat{{\mathbf {x}}}-\hat{{\mathbf {z}}}\hat{{\mathbf {z}}})$| for the zeroth order, while eq. (59) for the perturbation stresses become

where

given the relation eq. (43) for |$\mu _0$|⁠, and note that the dimensionless version of the basic state strain rate is

For the dimensionless equations, we drop the primes on perturbation quantities.

The dimensionless version of the governing equation for the stream function (eq. 65) becomes

The dimensionless version of the boundary conditions at the top |$z=0$| (eq. 67) are

The dimensionless version of the bottom boundary conditions at |$z=-1$|⁠, that is eq. (68), are

where

This system of equations, that is eq. (81) with boundary conditions (82) and (83), are solved with a standard pentadiagonal solver, an example of which is described in detail by (Bercovici et al. 2023, Appendix A1). The model is initiated by calculating the flow and stress fields associated with a uniform viscosity given by eq. (79) with eq. (76), and proceeding with the model iterations accordingly.

3.3 Model results

For this simple model, we explore only a few key effects on model behaviour, mostly investigating the influence of those parameters that are least well constrained, namely the grain damage numbers |$\mathcal {D}_i$| and |$\mathcal {D}_I$|⁠, as well as the phase mixing diffusion factor |$\chi$|⁠. We assume other parameters such as grain and interface coarsening |$\mathcal {C}_i$| and |$\mathcal {C}_I$|⁠, and drip buoyancy number |$\beta$| are better constrained and we fix those values to limit the range of parameter space (see Table 2). Moreover, we use values of damage and coarsening numbers extracted from the comparison to lab experiments of Bercovici et al. (2023), which provided important constraints on the relative size between the two damage numbers. However the experimental conditions that the model was compared to were at lower pressure, higher temperature and much higher strain-rates than typical of the lithosphere; thus the results should be interpreted accordingly. Again, the model is highly simplified and pedagogical and designed to give a working example of the upscaling framework. Moreover, note that the validity of the assumptions of small net lateral strain and linearized perturbation analysis are only accurate for |$t\ll 1$| and |$h\ll 1$|⁠, although we allow the solutions to evolve beyond these limits for the sake of illustration.

3.3.1 Drip evolution

One of the basic metrics of the mechanical response to the multiscale evolution of damage, mixing and weakening is the growth rate of the drip instability h (Fig. 5a). However, the development of anisotropic fabric from mesoscopic patch-scale mixing and microscopic grain-damage are also important, not only for establishing a potentially testable feature of lithospheric fabric, but also possibly a preferred direction of weakening that facilitates, for example, subduction initiation even before damage and weakening are ubiquitous (Bercovici & Mulyukova 2021, Figs 5b and c). We can gauge the evolution of the degree of anisotropy with a root-mean invariant of the Gary tensor |$\underline{{\boldsymbol{\Gamma }}}$| defined as

where |$\langle Q\rangle$| is the macroscale domain area average of Q (i.e. average over the domain |$x\in {-L/(2H),L/(2H)}$| and |$z\in {-1,0}$|⁠).

As noted in Section 2.4, the degree of anisotropy increases and decreases again as the mixing volume fraction |$\Phi$| goes from 0 to 1 (and, in principle, maximizes at |$\Phi ={1}/{_2}$|⁠). The magnitude of anisotropy also depends on the viscosity difference |$\Delta \mu = \mu _{\rm{\small {N}}}-\mu _{\rm{\small {S}}}$| (or non-dimensionally |$\Delta \nu$|⁠) between weak, mixed layers and strong, unmixed layers.

We examine how the lithospheric drip amplitude h and degree of anisotropy |$\bar{\Gamma }$| evolve through time versus grain-damage number |$\mathcal {D}_i$| and diffusive mixing coefficient |$\chi$|⁠, both of which are poorly constrained parameters; however, to limit the size of parameter space, we assume the grain and interface damage numbers are proportional according to |$\mathcal {D}_I=0.25\mathcal {D}_i$|⁠, as inferred by comparison to experiments (Bercovici et al. 2023). The healing numbers are held fixed since they are better constrained by laboratory measurements of grain-growth (e.g. Karato 1989), and specifically to |$\mathcal {C}_i=0.2$|⁠, |$\mathcal {C}_I=\mathcal {C}_i/2000$| (Bercovici et al. 2023, although note these arise from comparison to lab experiments at moderately low pressure, high stress and high temperature). The drip buoyancy number is held fixed to |$\beta =3$| (given typical mantle terrestrial numbers for density and lithospheric thickness and assuming a modest ridge push stress |$\tau _0$| of about 10 MPa). The lateral size of the domain is chosen such that |$L=2H$| and thus |$x\in [-1,1]$|⁠. The drip wavelength is |$\boldsymbol{k}=\pi$| (one cosine wave across the full horizontal domain). Otherwise we use canonical values of |$n=m=3$|⁠, |$p=2$|⁠, and |$q=4$| in eqs (27)–(29).

Varying damage number |$\mathcal {D}_i$| has multiple effects. First increasing |$\mathcal {D}_i$| causes a decrease in the background isotropic viscosity |$\mu _0$|⁠, or non-dimensionally |$\eta _0/2$|⁠, since greater damage reduces |$R_0$| and |$\eta _0 = R_0^m/(1+R_0^m)$|⁠. Secondly, increasing |$\mathcal {D}_i$| increases the viscosity contrast between mixed and unmixed zones (since Zener pinning amplifies grain-damage in mixed zones), thus leading to stronger influence on anisotropy. Thus the drip grows faster with increasing |$\mathcal {D}_i$| somewhat as expected given the influence on the mean viscosity |$\mu _0$| (or |$\eta _0/2$|⁠). The influence of viscous anisotropy on the drip evolution can be gauged by comparing the full anisotropic model to the isotropic case (wherein we enforce |$\nu =\Gamma _{ij}=0$|⁠) which is equivalent to single-phase damage occurring without any mixing at the mesoscale (see Fig. 5a, black curves). The difference between drip growth rates for the anisotropic and isotropic cases is significant but not remarkable for the short times of the calculations; however the curves diverge and at much larger times than is reasonable for the linear approximation, they will potentially diverge further and allow significantly faster drips for the anisotropic media.

More remarkable than drip growth is that the maximum degree of anisotropy |$\bar{\Gamma }$| increases several orders of magnitude for only a 1 or 2 order of magnitude increase in |$\mathcal {D}_i$|⁠, which is due to the large grain size sensitivity of viscosity (both in the mean viscosity |$\mu _0$| – which decreases with increasing |$\mathcal {D}_i$| – and viscosity contrast |$\Delta \mu$| – which increases with |$\mathcal {D}_i$| – and given that both affect |$\Delta \nu = \Delta \mu /\mu _0$| in the relation for |$\underline{{\boldsymbol{\Gamma }}} = \Delta \nu \underline{\mathbf {J}}$|⁠). But only at larger |$\mathcal {D}_i\approx 10$| is the anisotropy significant enough to influence the net viscosity.

Anisotropy, specifically |$\bar{\Gamma }$|⁠, is also influenced by the mixing diffusion coefficient |$\chi$| (Fig. 5c). Mixing is generally very slow anyway, and we expect |$\chi\ll 1$|⁠; but given the stretching of the phase field during deformation, diffusive mixing eventually dominates and eliminates heterogeneity, causing both |$\phi ^{\prime }\rightarrow 0$| (or |$\phi \rightarrow \phi _0$|⁠) and for |$\Phi \rightarrow 1$| (complete ubiquitous mixing), which corresponds to isotropy. For larger |$\chi$|⁠, diffusive mixing is reached more rapidly; thus the interval for which |$\bar{\Gamma }$| goes from zero (i.e. an isotropic state before mixing ensues and |$\Phi =0$|⁠), to a maximum value with peak anisotropy (near |$\Phi ={1}/{_2}$|⁠), and back to zero again (as |$\Phi \rightarrow 1$|⁠) occurs sooner, lasts more briefly, and has smaller peak value of |$\bar{\Gamma }$|⁠. As |$\chi$| becomes much smaller, anisotropy starts later, but lasts longer and reaches higher peaks, that have a stronger influence on the viscosity; with smaller |$\chi$|⁠, the phase heterogeneity becomes stretched even more into an anisotropic laminated structure before weakening.

3.3.2 Flow and viscosity fields evolution

The evolution of various fields at the macroscale are informative for inferring the influence of the meso- and microscales (Figs 6, 7 and 8). The macroscopic fields that are diagnostic of system evolution are stream function |$\psi$|⁠; mean grain size |$\bar{R}$|⁠; the unique components of the stress, Gary and viscosity tensors |${\underline{{\boldsymbol{\tau }}}}$|⁠, |$\underline{{\boldsymbol{\Gamma }}}$| and |$\underline{\mathbf {M}}$|⁠; the drip depth h; and a spatial sampling of the three phase field vectors |${\mathbf {k}}_j$|⁠, the initial values of which are in eq. (7) (for |$N=3$|⁠) but evolve according to eq. (15). We again vary |$\mathcal {D}_i$| and |$\chi$| (as in Fig. 5).

The shape of the stream function |$\psi$| changes little with varying control parameters, but that is largely because we assume it has a single-mode sinusoidal structure; thus only the amplitude of |$\psi$| increases as |$\mathcal {D}_i$| increases and the isotropic viscosity decreases.

The mean grain size |$\bar{R}$| is reduced where the vertical tensile stress (⁠|$\tau _{zz}\equiv -\tau _{xx}$|⁠) of the drip augments the vertical tensile background stress (imposed by the lateral ‘ridge push’ compression) and leads to a vertical tongue of small grain size and low viscosity right above the nascent drip or slab, which potentially augments subduction initiation. The effect is more pronounced with increasing damage number |$\mathcal {D}_i$|⁠, as expected (see also Paczkowski et al. 2012).

The stress fields |${\underline{{\boldsymbol{\tau }}}}$| also show the enhanced tensile stress from the drip, but lower shear stresses due to return flow following the drip. The viscosity |$\underline{\mathbf {M}}$| is lower where grain sizes are small, in a vertical tongue over the drip. The off-diagonal components of the viscosity field |$M_{xz}$| are always smaller than the diagonal ones |$M_{xx}$|⁠. The Gary tensor |$\underline{{\boldsymbol{\Gamma }}}$| is negligibly small for low |$\mathcal {D}_i=0.2$|⁠, and is barely reflected in the viscosity field. But as |$\mathcal {D}_i$| increases, the Gary tensor components get much larger (also seen in Fig. 5b) and are more reflected in the viscosity field. Specifically, while |$M_{xz}$| always mirrors |$\Gamma _{xz}$|⁠, the central vertical viscosity minimum in |$M_{xx}$| becomes narrower and more closely mirrors the extrema in |$\Gamma _{xx}$| as |$\mathcal {D}_i$| increases from 0.2 to 10 (Figs 6–8). The three phase field wave vectors |${\mathbf {k}}_j$| start off with the hexagonal arrangement (each |$60^\circ$| from each other) and then get deflected in the compressive stress direction |$\pm \hat{{\mathbf {x}}}$|⁠, meaning the phase heterogeneity aligns in vertical lenses (for example as shown in Fig. 3a); but as the drip grows the orientations become a little more scattered especially near the base where the most compressive directions are deflected by the drip surface motion.

Variation in the diffusion coefficient |$\chi$| has little bearing on the drip evolution, but has notable influence on the development of anisotropy, as also illustrated in Fig. 5c). Specifically, the peak anisotropy reflected in |$\underline{{\boldsymbol{\Gamma }}}$| occurs sooner and is of lower magnitude for increased |$\chi$|⁠, and the diagnostic fields at the times of these peaks are markedly distinct for different |$\chi$| (Fig. 9). Because larger |$\chi$| causes faster diffusive mixing, the mixing zones develop more rapidly, leading to a peak anisotropy arriving sooner, but then total mixing and the obliteration of anisotropy happens more rapidly too. For very small values of |$\chi$|⁠, mixing zones take longer to develop, but then take much longer to be erased, and in the meantime the stretching of the mixed zones amplifies the level of anisotropy. For the same reasons, cases with large |$\chi$| show little effect of anisotropy in the net viscosity field |$\underline{\mathbf {M}}$|⁠, while cases with small |$\chi$| show a stronger influence of anisotropy in the viscosity.

4 DISCUSSION AND CONCLUSION

The primary goal of this paper is to provide a practical framework for solving and implementing the multiscale physics of coupled mineral phase mixing and grain-damage into geodynamic models. Specifically, deformation driven mixing between mineral phases via grain migration of one phase along the opposite phase’s grain-boundaries (either mechanically or chemically) occurs at the small (centimetre) scale of petrological heterogeneity. The enhancement of grain-damage and development of mylonitic weak zones occurs within the mixed zones between phases, but this involves microscale grain size processes, that is within and between grains, such as dislocation movement and accumulation and dynamic recrystallization, grain-boundary migration and grain-growth, and Zener pinning between grains of different minerals. All these centimetre to micron scale processes influence the strength and mechanical behaviour of the mantle and lithosphere at the kilometre to hundred kilometre scale. Capturing the physics at these small meso (heterogeneity) and micro (grain) scales to inform flow and deformation at the macro (tectonic) scale is both a challenge and a necessity. Clearly, if each point on a macroscopic space is itself a finite domain or a ‘mesoscopic patch’ that also involves multidimensional transport such as diffusive mixing (and within each patch, the mixed and unmixed zones have different grain evolution), then solving the physics at all scales perfectly becomes untenable.

We therefore propose a framework to model the small meso- and microscale processes with efficient and rapid solutions that capture the essential physics with little computational cost. Specifically, the potentially time-consuming grain-mixing diffusion at the mesoscale patch is resolved into a simple analytic solution by tracking a minimum number of cosine modes in the phase field that is sufficient to model heterogeneity, and transforming into an appropriate coordinate system in the patch that yields an analytic solution for how heterogeneity evolves during simultaneous diffusion and deformation. This solution provides information about the structure and volume fraction of mixing zones. The microscale grain-scale processes are then also solved quickly by attributing two-phase grain-damage with pinning within the mixing zones, and single phase damage in the unmixed zones. These grain-damage calculations result in grain sizes and viscosity estimates in the mixed and unmized zones, which are then used to infer a bilaminate viscosity structure (wherein the mixed zones readily develop into weak sheets of small grain sizes—i.e. mylonitic layers) that in the macroscopic reference frame manifests as an anisotropic viscosity field. This final viscosity field is then communicated to the geodynamic model to infer macroscopic flow and stress, which can then be imposed again on the meso- and microscopic processes. Each scale obeys the same time and stress scale (i.e. imposed by forcing and deformation) but each has a unique length scale. Non-dimensionalizing by the appropriate length scale also leads to easily resolved solutions at the various scales.

The application of this framework to a simple example of a lithospheric layer under ridge-push-like lateral compressive stress, and a Rayleigh–Taylor like drip instability, demonstrates how the framework can be implemented, but also the effect of the subscale processes on the development of viscous anisotropy and subduction initiation. The deforming model lithosphere experiences both weakening in the isotropic viscosity field (i.e. the average of all the microscale weakening), and also development of anisotropy due to the mesoscale mixing with microscale grain damage. For the small-amplitude drip calculations, anisotropy makes only a modest contribution to the drip instability, although it potentially becomes much larger for finite-amplitude drips. Weakening and anisotropic fabric in a tongue above the drip imply that the lithosphere above a drip or a nascent subduction could be weakened to promote subduction initiation, although it is premature to draw too many conclusions from a linear, small-amplitude model.

The drip rate is, as expected, accelerated with increasing grain-damage. But the degree of anisotropy is also very sensitive to grain-damage, increasing orders of magnitude with only a moderate increase in grain-damage. But since anisotropy in the model is due to the alignment and weakening of mixing zones, it is impermanent so long as deformation and mixing persists. Specifically, before mixing ensues there are no mixed zones and the medium is isotropic and strong. When mixing has persisted so long that the medium is uniformly mixed, then again there are no mixed bands and the medium is isotropic again but weak. Only at intermediate times does the medium reach peak viscous anisotropy from these processes; however, given that mixing by grain migration is intrinsically slow and limited by molecular diffusion (either chemical diffusion or diffusion creep) then the longevity of anisotropy might be extensive even if it is ephemeral.

This paper proposes a framework for an efficient solution to coupled, multiscale physical processes, but much work remains to be done to explore the effects to this type of model and approach. The efficiency of solving the mesoscale diffusion system relies on an analytical solution assuming constant diffusivity; allowing for a variable diffusion coefficient |$\mathcal {K}$| that is possibly a function of grain size and phase fraction (Bercovici & Mulyukova 2018, 2021) would require numerical solutions, but can be done in such a way as to keep them as simple ordinary differential equations (no more difficult to solve than the grain size evolution laws). The analytical solution also depends on the velocity being linear across the patch and thus of homogeneous strain, even though there are weak and strong bands; this means the weakening is due to stress reduction, not strain-rate enhancement, which constrains the length scale of strain localization to be between, not within, patches. The assumption of bilaminate weakening is also worth relaxing to allow, for example, elliptical (in 2-D) or ellipsoidal (in 3-D) weak inclusions; however given the straining of weak layers, laminar structures are probably not a bad assumption at least at large strains.

Extending the upscaling scheme to three-dimensions is also an obvious next step and is outlined for completeness in Appendix  A. Moreover, the framework can be extended to cases where phases have different viscosities and anisotropy is due to stirring of phases and reaching stress percolation limits (localization within the weak phase), even without diffusive mixing and mylonitization in mixed zones. The interaction of implicitly weak phases with phase mixing is also an area in need of exploration both theoretically and experimentally (e.g. Cross et al. 2020). Finally, a fully non-linear model application of the framework to a lithospheric drip would be desirable, as would even general convective evolution of a cold thermal boundary layer, to name a few examples. In the end, the potential to explore the coupling of these multiscale processes hopefully leads to new insight in the collaboration between mineral physics, rock deformation and geodynamic modelling of the mantle and lithosphere.

Acknowledgement

Many thanks to the editor Henri Samuel, and to Juliane Dannberg and Laurent Montesi for their comprehensive and thoughtful reviews.

DATA/CODE AVAILABILITY

All original results are numerical. Matlab^® codes for the numerical models and supporting routines and sample data for Figs 5–9 are available on Zenodo (doi:10.5281/zenodo.11479173).

References

APPENDIX A: NOTES ON THE 3-D FRAMEWORK

The framework for upscaling presented in the main text assumes plane-stress and 2-D flow, which is reasonable for applications such as subduction initiation. The framework is extendable to 3-D, which is less elegant and more cumbersome than in 2-D, but is nevertheless tractable. For the sake of completeness, we outline the necessary steps to pose the upscaling framework in three dimensions.

The primary difference between the 2-D and 3-D formulation is in the mesoscale level of the model, i.e. in terms of describing diffusive grain-mixing in the mesoscale patch (or volume), the fabric of mixed zones, and the resultant anisotropic viscosity. The 3-D macroscopic formulation is also different although only by dint of requiring an extra degree of freedom in the flow field; although the approach is standard we will outline the macroscopic formulation in 3-D as well. The microscopic grain-scale formalism does not change from the 2-D framework, other than using the 3-D expressions for the stress and strain-rate invariants as well as the mixing work in the various grain-boundary and interface damage terms.

A1 Principal stress directions and compressive components

Given a 3-D stress tensor |${\underline{{\boldsymbol{\tau }}}}$|⁠, the principal stresses are the solutions to the cubic polynomial

where |$\tau ^2 ={1}/{_2}{\underline{{\boldsymbol{\tau }}}}:{\underline{{\boldsymbol{\tau }}}}$| and |$\mathrm{det}({\underline{{\boldsymbol{\tau }}}})$| are the second and third invariants of |${\underline{{\boldsymbol{\tau }}}}$| (and, having assumed incompressibility, the first invariant |$\mathrm{trace}({\underline{{\boldsymbol{\tau }}}}) = \tau _{kk} = 0$|⁠). Eq. (A1) is already in the form of a depressed cubic polynomial that has 3 real analytic roots for s, and these are the three principal stresses |$\tau _1, \tau _2$| and |$\tau _3$| given by

(see Selby 1973). The principal stresses are ordered such that |$\tau _1$| is the most positive and tensile, and |$\tau _3$| is the most negative and compressive; |$\tau _2$| can be either positive or negative. (We also recognize that the subscripted stresses |$\tau _k$| are similar to |$\tau _i$| for stresses in phases |$i=1$| or 2 discussed in the main text, specifically in Section 2.3 about processes at the micro-scale grain size space. However, in this appendix we do not refer to stresses in the phases, and rather than introduce yet more symbols, we retain the |$\tau _k$| for the principal stresses and hope it does not cause excessive confusion or convulsions.)

The three orthonormal eigenvectors (principal directions) associated with |$\tau _1, \tau _2$| and |$\tau _3$| are |$\hat{{\mathbf {l}}}, \hat{{\mathbf {m}}}$|⁠, and |$\hat{{\mathbf {n}}}$|⁠, respectively. (Note that |$\hat{{\mathbf {n}}}$| is the same as in the 2-D formulation, that is in the principal compressive stress direction, but |$\hat{{\mathbf {m}}}$| is not, since now |$\hat{{\mathbf {l}}}$| is the principal tensile direction.) In general, these eigenvectors have the form of

The compressive component of |${\underline{{\boldsymbol{\tau }}}}$| always involves |$\tau _3\hat{{\mathbf {n}}}\hat{{\mathbf {n}}}$| and may or may not involve the less compressive component |$\tau _2\hat{{\mathbf {m}}}\hat{{\mathbf {m}}}$|⁠, depending on the sign of |$\tau _2$|⁠. We can write the compressive stress as

where we define |$\tilde{\tau }_2 = {1}/{_2}(1-\mathrm{sign}(\tau _2))\tau _2$|⁠.

A2 Mixing diffusion

To model grain-mixing as a diffusive process in a mesoscale patch, we again do a Galilean transformation to the frame of reference of the mean velocity of the patch, and then rotate into the principal stress directions described in (A3). As in the 2-D model, the mixing diffusion equation is

where |$\tilde{{\mathbf {v}}}$| is the transformed velocity, relative to the patch mean (or equivalently central) velocity and rotated to the stress principal directions. Taking the expression for |${\underline{{\boldsymbol{\tau }}}}_c$| from (A4), (A5) becomes

As in the 2-D formulation, we assume that, in the principal stress coordinate system, the deformation is described by pure shear and, in which the strain-rate tensor expressed in the principal stress frame is

where |$\dot{e}_1+\dot{e}_2+\dot{e}_3=0$|⁠, and |$\dot{e}_1, \dot{e}_2$| and |$\dot{e}_3$| are the principal strain-rates of the macroscopic strain-rate |$\dot{\underline{\mathbf {e}}}$|⁠. (See also Section 2.2.3 for discussion of this assumption.) The vorticity tensor is

where the three |$\Omega _j$| are the independent components of |$\underline{{\boldsymbol{\Omega }}}$| in the the principal stress coordinate system (i.e. determined from the macroscopic velocity field in the |$(x,y,z)$| reference system and transformed into the principal stress system.). Note that |$\underline{{\boldsymbol{\Omega }}}$| obeys the same skew symmetry in any coordinate system. The |$\Omega _j$| are also proportional (within a factor of |$\pm {1}/{_2}$|⁠) to the components of the negative of vorticity in the |$\hat{{\mathbf {l}}}$|⁠, |$\hat{{\mathbf {m}}}$| and |$\hat{{\mathbf {n}}}$| directions, respectively. The transformed, patch-centred velocity field is assumed to obey a linearized expression as in eq. (10), which leads to

The representation of |$\phi$| in cosine structures is similar to that in eq. (6) except the position vector |$\tilde{{\mathbf {x}}}$| (relative to the centre of the patch), and the wave vector |${\mathbf {k}}$| are 3-D, the latter according to |${\mathbf {k}} = k_l\hat{{\mathbf {l}}}+k_m\hat{{\mathbf {m}}}+k_n\hat{{\mathbf {n}}}$|⁠. Specifically, instead of eq. (6), we can write

and before shearing and diffusion take place, all the |$\breve{\phi }_{ij}$| are the same, such that for example |$\breve{\phi }_{ij}=\frac{1}{N(N+1)/2}$|⁠, and the wave vectors obey

where |$k=2\pi /\lambda$|⁠, and again |$\lambda$| is the characteristic wavelength between phase anomalies. For example, if |$N=2$| there are |$N(N+1)/2=3$| orthogonal |${\mathbf {k}}_{ij}$| that describe a cubic structure. The above expression for |${\mathbf {k}}_{ij}$| is somewhat arbitrary and is shown as an example; other expressions can be used, although at the mesoscopic scale of |$\delta=\mathsf {n}\lambda$|⁠, the heterogeneity should be coarse (carrying few modes) and quasi-periodic.

The evolution and distortion of the wave vector |${\mathbf {k}}$| (dropping the |$ij$| subscripts for now) still obeys eq. (13), which, given eq. (A9), leads to the system of coupled equations

The general solution of eq. (A12) is

where |$\imath =l,m$| or n, and the three eigenvalues |$\alpha _{\jmath }$| are the three roots to the cubic polynomial characteristic equation arising from

The nine constants |$A_{\imath \jmath }$| are determined from the three initial conditions that |$k_{\imath }(t_0)=\sum _{\jmath }A_{\imath \jmath }=K_{\imath }$|⁠, for example as in eq. (A11), and that two of the equations in eq. (A12) must be obeyed for each of the three |$\alpha _{\jmath }$| [the third remaining relation in eq. (A12) is already satisfied by the eigenvalue relation (eq. A14) for |$\alpha _{\jmath }$|], leading to six more constraints.

Assuming that |$\mathcal {K}$| is uniform over the patch (as is the macroscopic stress), the diffusion equation for each Fourier mode amplitude becomes

for which the analytic solution can be readily obtained using (A13). The full phase field within the patch still then evolves according to eq. (21) but summing over i and j as in eq. (A10), and with |${\mathbf {k}}$| and |$\tilde{{\mathbf {x}}}$| in 3-D.

A3 Fabric tensor

Given an analytic solution for |$\phi$| in each patch, the fabric tensor derived from the mixing flux is

But we can define the following components of the tensor for convenience

The two eigenvalues of the fabric tensor are

where the root |$F_+$| is associated with the largest mixing flux (and in the case where there is only one compressive stress direction and |$\tilde{\tau }_2=0$|⁠, then |$F_+= F_{nn}$| and |$F_-=0$|⁠). The eigenvectors associated with |$F_+$| and |$F_-$| are

respectively. Since the vector |$\hat{{\mathbf {q}}}$| is associated with the fastest mixing flux, it is normal to the mixing bands, while |$\hat{{\mathbf {p}}}$| is parallel to them (again consider the case when |$\tilde{\tau }_2=0$|⁠), and these will be used next to construct the anisotropic viscosity field. Moreover, note that

is the mixing work term to be used in the grain-damage relations [i.e. use |$\mathcal {F}$| in lieu of |$\mathcal {F}_n$| in eq. (27b)].

Finally, we note that the mixing zone volume fraction |$\Phi$| arises from the same formalism as in the 2-D model, specifically eqs (25) and (26), except for summing over i and j

A4 Anisotropic viscosity

For simplicity, and consistency with the 2-D model, we assume the mixing bands effectively form bilaminate weak structures, although more complex geometries can also be invoked (such as ellipsoidal bands; see Treagus 2002, 2003; Fletcher 2004; Montési 2007; Dabrowski et al. 2012).

The grain sizes and viscosities in the mixed and unmixed zones follow the same two-phase grain-damage relations as in the 2-D model (see Section 2.3), and the effective viscosities |$\mu _{\rm{\small {N}}}$| and |$\mu _{\rm{\small {S}}}$| for deformation normal and in shear parallel to the laminate weak layers follows the same relation as in eq. (34).

As derived above, the unit vectors normal and parallel to the laminate weak layers are |$\hat{{\mathbf {q}}}$| and |$\hat{{\mathbf {p}}}$|⁠, respectively, and both are in the same plane as the principal stress directions |$\hat{{\mathbf {n}}}$| and |$\hat{{\mathbf {m}}}$|⁠, and normal to |$\hat{{\mathbf {l}}}$|⁠. The five independent stress components oriented with respect to the laminate layers are

where |$\tau _{ll} + \tau _{qq} + \tau _{pp} = 0$|⁠. The transverse isotropic constitutive law in 3-D is

where the strain-rate components have an analogous relation to eq. (A21). We assume that any shear stress involving a traction parallel to and on the surface of a weak laminate layer involves the weaker viscosity |$\mu _{\rm{\small {S}}}$|⁠, while all other stresses involve the stronger viscosity |$\mu _{\rm{\small {N}}}$|⁠.

The stress components (eq. A21) (and implicitly the strain-rate components) are related to the original Cartesian components through the relation between |$\hat{{\mathbf {q}}}$| and |$\hat{{\mathbf {p}}}$| to |$\hat{{\mathbf {n}}}$| and |$\hat{{\mathbf {m}}}$| in eq. (A19) and then the relation between |$\hat{{\mathbf {l}}}$|⁠, |$\hat{{\mathbf {m}}}$| and |$\hat{{\mathbf {n}}}$| to |$\hat{{\mathbf {x}}}, \hat{{\mathbf {y}}}$| and |$\hat{{\mathbf {z}}}$| in eq. (A3). For example, |$\tau _{qq}$| is linearly related to all five independent Cartesian stress components through

and similarly for the other stress (and strain-rate) components, leading to the overall linear relation

where |$\tau _{xx}+\tau _{yy}+\tau _{zz}=0$|⁠, and similarly for the strain-rate tensor. The transverse isotropic constitutive law becomes, as in the 2-D model,

where the inverse of the matrix |$[Q]$| is presumably the transpose of |$[Q]$| given that it is built on Cartesian transformations, although we defer proof of that to the young and ambitious (and patient). This final relation is, as with the 2-D model, written simply as

where |$\underline{\mathbf {M}}$| is the anisotropic viscosity field given by the product of 5×5 matrices in eq. (A25).

A5 Macroscopic flow

The governing equations for macroscopic flow is no different from that in the 2-D model, that is (47)–(49). Given incompressibility, there are only two solenoidal flow potentials given by

which are the toroidal and poloidal fields, respectively, and that also prescribe the strain-rate field. Thus only two components of the force equation are necessary to solve for these potentials. The first arises from |$\hat{{\mathbf {z}}}\cdot \boldsymbol{\nabla }\times$| the force-balance eq. (48), which in terms of stress components leads to

The second relation arises from |$\hat{{\mathbf {z}}}\cdot \boldsymbol{\nabla }\times \boldsymbol{\nabla }\times$| of eq. (48) and leads to

Given the anisotropic constitutive laws (eq. A26), these equations provide relations for the toroidal and poloidal potentials |$\Psi$| and W, as also illustrated in the 2-D model for the solenoidal stream function |$\psi$|⁠. It is straight-forward to show that in the isotropic limit, (eq. A28) and (eq. A29) reduce to

where |$\nabla _h^2 = \frac{\partial ^2 }{\partial x^2}+\frac{\partial ^2 }{\partial y^2}$|⁠. As with the 2-D framework, the calculated macroscopic flow field is used to determine the stress field from eq. (A26), which then dictates the new principal stresses and stress directions, with which to update all the processes at the mesoscopic (patch) and microscopic (grain) scales.

APPENDIX B: TWO-DIMENSIONAL DRIP MODEL BOUNDARY CONDITIONS

The model lithospheric layer in the 2-D example provided in the text (Section 3) is assumed to be much more viscous than the underlying asthenosphere, and certainly more so than the ocean or atmosphere above; thus the layer is bounded by relatively inviscid media such that its top and bottom boundaries can be treated as free-slip (shear-stress free). Moreover, we assume that the density contrast between the layer and the overlying medium is very large and that deflections of that boundary are negligible; thus the top is also treated as impermeable (zero vertical velocity). The bottom boundary however supports deflections of thickness |$h(x,t)$|⁠, which influence the stress and velocity boundary conditions.

B1 Top boundary

At the top boundary, where |$z=0$|⁠, the impermeable condition requires that the vertical velocity is zero. Since the zeroth order vertical velocity, due to lateral compression, is |$w_0 = -\dot{e}_0z$|⁠, it is already zero at |$z=0$|⁠. Thus the perturbation velocity is likewise zero:

The free-slip condition requires zero shear stress on the top boundary, which being approximately flat and undeformed, and using eqs (59b) and (62), leads to

B2 Bottom boundary

At the bottom boundary, where |$z=-H-h$|⁠, the zeroth order vertical velocity is |$w_0 = -\dot{e}_0H$|⁠, which implies that H itself is time dependent and goes as |$H = H_0(1 + \dot{e}_0t)$|⁠. Thus, strictly speaking, the model domain itself deforms, shrinking laterally and expanding vertically. However, for simplicity (and trying to adhere to a simple demonstration), we only assume |$\dot{e}_0 t\ll 1$|⁠, so that we can approximated H as constant; should one consider times |$t\gt 1/\dot{e}_0$|⁠, the spatial domain should be made time dependent.

The perturbation in vertical velocity is

(Strictly speaking, this boundary condition should involve the evaluation of |$\breve{\psi }$| at |$z=-(H+h)$|⁠, but expanding |$\breve{\psi }(-H-h)$| around |$z=-H$| only adds a second order perturbation that we would drop anyway.) However, this boundary condition provides an evolution equation for h, and not necessarily a boundary condition for the flow field. Given that the bottom boundary is deformable, we need to consider continuity of all stress components across it. The traction vector (force per area vector) across the bottom boundary is |$(-P\underline{\mathbf {I}}+{\underline{{\boldsymbol{\tau }}}})\cdot {\hat{{\boldsymbol{h}}}}_\perp$|⁠, accounting for the full stress on the boundary, and where |${\hat{{\boldsymbol{h}}}}_\perp$| is the unit normal to the surface given by |$z+H+h=0$|⁠, such that

(see Fig. 4). There are zeroth and first order terms in pressure and stress such that |$P = P_\text{h}(z) + P^{\prime }({\mathbf {x}},t)$|⁠, in which |$P_\text{h}$| is the hydrostatic pressure and |$P^{\prime }$| the perturbation due to viscous motion, and |${\underline{{\boldsymbol{\tau }}}}$| is given by eq. (57); in this case the traction vector on the boundary is

where we drop the second order products of perturbations in stress and boundary deflection. The normal stress at the boundary is

(where again recall |$\tau ^{\prime }_{xx}=-\tau ^{\prime }_{zz}$|⁠). The shear stress at the boundary is the component of the traction vector that is parallel to the boundary, or in the direction of the tangent unit vector |${\hat{{\boldsymbol{h}}}}_{||}= -\hat{{\mathbf {z}}}\frac{\partial h}{\partial x} + \hat{{\mathbf {x}}}$| (wherein |${\hat{{\boldsymbol{h}}}}_\perp \cdot {\hat{{\boldsymbol{h}}}}_{||}=0$|⁠; see Fig. 4), and is given by

The continuity of shear stress across the boundary requires a jump condition |$\left[\tau ^{\prime }_{xz}-2\tau _0\frac{\partial h}{\partial x}\right]_-^+ = 0$|⁠, where |$[..]_-^+$| denotes the difference between a quantity evaluated just above (⁠|$+$|⁠) and below (⁠|$-$|⁠) the boundary. Since the medium below the boundary is inviscid it supports no deviatoric stresses and thus at |$z=-(H+h)$| we require |$\tau ^{\prime }_{xz} = 2\tau _0\frac{\partial h}{\partial x}$|⁠. We project this condition to the position |$z=-H$| by considering perturbations to it by the boundary deflection; however, this only leads to additional second order terms (e.g. by expanding |$\tau ^{\prime }_{xz}(-H-h)$| to first order in h leads to |$\tau ^{\prime }_{xz}(-H)$| plus second order terms). Thus to first order, the zero shear-stress, or free-slip, condition on the bottom boundary is

Physically this states that the perturbation shear stress on the boundary must balance the effect of the zeroth order normal stresses pushing against the slope of the boundary.

The continuity of normal stress across the boundary requires the jump condition |$\left[-P_\text{h}+\tau _0 - P^{\prime } + \tau _{zz}^{\prime }\right]_-^+ = 0$| at |$z=-(H+h)$| which leads to

where |$\Pi _\text{h}$| is the hydrostatic pressure in the inviscid underlying medium, where all other quantities like deviatoric stress and viscous pressure perturbations are zero. The above relation can be expanded to first order around |$z=-H$| to yield

where we have exploited the hydrostatic balance |$\frac{\partial P_\text{h}}{\partial z}=-\rho g$|⁠, and recall that the underlying medium is |$\Delta \rho$| less dense than the lithospheric layer with density |$\rho$|⁠. The zeroth and first order conditions are satisfied separately; thus |$P_\text{h}-\Pi _\text{h} = \tau _0$| allows for a jump in hydrostatic pressure due to the lateral compression of the layer, and the perturbation in normal stress satisfies

These stress conditions are, however, only useful when put in terms of the stream-function |$\psi$|⁠. The free-slip condition (eq. B8) along with eqs (59b) and (62), leads to

where we have used the definition |$\dot{e}_0 = \tau _0/(2\mu _0)$|⁠, and recall |$h = \breve{h}\cos (\boldsymbol{k}x)$|⁠, while |$(\psi , \Gamma _{xz}) = (\breve{\psi }, \breve{\Gamma }_{xz})\sin (\boldsymbol{k}x)$|⁠.

To put the normal stress condition in terms of the stream function |$\psi ^{\prime }$| only, we take |$\frac{\partial }{\partial x}$| of eq. (B11) and use the x component of the force balance eq. (48) to eliminate |$P^{\prime }$|⁠, which yields

where we used |$\tau ^{\prime }_{zz}=-\tau ^{\prime }_{xx}$|⁠. Using (59) with eq. (62), the above relation becomes

which can be reorganized to

The above equation can be reposed in terms of the sinusoidal modal amplitudes

and again recall |$(h,\nu ,\Gamma _{xx}) = (\breve{h}, \breve{\nu },\breve{\Gamma }_{xx})\cos (\boldsymbol{k}x)$|⁠, while |$(\psi , \Gamma _{xz}) = (\breve{\psi }, \breve{\Gamma }_{xz})\sin (\boldsymbol{k}x)$|⁠.