# Two-phase damage and plate generation in a 2-D model of mantle convection

William Landuyt, David Bercovici, Yanick Ricard

► To cite this version:

William Landuyt, David Bercovici, Yanick Ricard. Two-phase damage and plate generation in a 2-D model of mantle convection. Geophysical Journal International, 2008, 10.1111/j.1365-246X.2008.03844. hal-00361567

HAL Id: hal-00361567

https://hal.science/hal-00361567v1

Submitted on 17 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Plate generation and two-phase damage theory in a model of mantle convection

W. Landuyt,<fnref n="1" /> D. Bercovici<fnref n="1" /> and Y. Ricard<fnref n="2" />

Accepted 2008 May 6. Received 2008 February 2; in original form 2007 June 1

## SUMMARY

The formation of narrow, rapidly deforming plate boundaries separating strong plate interiors are integral components of the generation of plate tectonics from mantle convection. The development of narrow plate boundaries requires the interaction of a non-linear rheology and convection. One such non-linear rheology is two-phase damage theory which employs a non-equilibrium relation between interfacial surface energy, pressure and viscous deformation, thereby forming a theoretical model for void generation. Two-phase damage theory was recently extended to allow for deformational work to increase the fineness (reduce the grain size) of the matrix phase. We present results testing two-phase damage theory in a 2-D convectively driven system where we allow for (1) pure void-generating damage, (2) pure fineness-generating damage and (3) combined void- and fineness-generating damage. Pure void-generating damage is found to be unsuccessful at producing plate-like features. Fineness-generating damage is successful at inducing plate-like behaviour in certain circumstances, including increasing viscosity sensitivity to fineness and certain regimes of damage input and healing rate. Cases with combined void- and fineness-generating damage produce significantly more localization than the end-members due to the apparent increase of deformational work input into fineness generation. The interaction of microcracks and grain size reduction in two-phase damage theory suggests a rheological model for shear localization necessary for the formation of plate tectonic boundaries.

**Key words:** Planetary tectonics; Dynamics of lithosphere and mantle.

## 1 INTRODUCTION

Plate tectonics is the unifying theory of geology, yet much remains to be understood about the development of plate tectonics on Earth as well as its absence on the other terrestrial planets (Bercovici et al. 2000). Plate tectonics on Earth has been determined to have been in operation for at least two billion years and may well have been in operation much earlier (Cawood et al. 2006). These observations place an important constraint on numerical models in that Earth-like simulations of plate motion must be non-sporadic and long-lasting. The collection of studies labelled as plate generation refer to the attempt to understand how plate tectonics (or the surface manifestation of mantle convection) arise self-consistently within planetary convection, as opposed to viewing plate tectonics and planetary convection as separate entities. One of the most prominent features of plate tectonics that plate generation studies attempt to address is the existence of broad, strong plate interiors with weak, rapidly deforming boundaries (Weinstein & Olson 1992; Moresi & Solomatov 1998). The rapidly deforming boundaries tend to be very narrow features and are synonymous with shear localization. Also, the existence of strike-slip boundaries and toroidal motion remain

some of the most challenging features for plate generation studies to explain (Bercovici 2003). Initiation of subduction (King 2001) and the creation of passive ridges (Ricard & Froidevaux 1986; Tackley 2000b) are also vital aspects of plate tectonics to be understood.

The features of plate tectonics discussed above imply that the plates do not deform according to a simple Newtonian rheology. Experimental work on upper mantle rocks has shown that a significant portion of the upper mantle deforms according to dislocation creep with a power law of approximately three (e.g. Karato & Wu 1993). Incorporation of non-Newtonian rheologies into plate generation studies has not shown a significant improvement in generating plate-like behaviours, especially in producing toroidal motion (Christensen & Harder 1991; Weinstein & Olson 1992; Bercovici 1993, 1995). Implementation of more exotic rheologies (e.g. stick-slip rheologies) has had more success at generating plate-like behaviour (Bercovici 1993; Tackley 1998, 2000b), but the relevance of such rheologies to the actual lithosphere is questionable (Bercovici 2003). The longevity of plate boundaries (and their ability to be reactivated) implies that the low viscosity zone that constitutes the boundary is a time-dependent quantity that only decays after very long times (Gurnis et al. 2000). Therefore, tracking the time

<footnote n="1">¹Department of Geology and Geophysics, Yale University, New Haven, CT, USA. E-mail: william.landuyt@yale.edu</footnote>

<footnote n="2">²Laboratoire des Sciences de la Terre, CNRS, Ecole Normale Superieure de Lyon, Lyon, France</footnote>
evolution of the viscosity (or its controlling parameter) is imperative for plate generation studies. Previous studies have attempted to incorporate a time-dependent damage parameter (e.g. void or defect creation) to induce severe weakening in the plate that would address some of the deficiencies of the previous models (Bercovici 1998; Tackley 2000b; Auth *et al.* 2003; Ogawa 2003). These models have been able to generate narrow shear zones, but the connection between the assumed damage parameter and the underlying physics of damage creation is often prescribed rather than derived (Bercovici 2003).

Two-phase dynamics has been richly studied in geophysics, with many people considering its application to problems in melt-dynamics (McKenzie 1984; Spiegelman 1993a,b; Katz *et al.* 2006). Two-phase damage theory was developed in order to take a first-principles approach to studying the partitioning of deformational work into both dissipative heating and the generation of surface energy (Bercovici *et al.* 2001a). The generation of surface energy represents internal damage within the medium, and the manifestation of surface energy in the medium is through the creation and/or expansion/dilation of voids. Two-phase theory is, therefore, developed in order to track both the rock and void phases. Since voids are suggested to represent damage or weakness in the material, tracking the voids essentially boils down to tracking the damage in the system. The theory was generalized to allow for deformational work to go into generating fineness (grain size reduction) in the rock phase (or essentially pulverizing the rock) in addition to void generation, and pulverizing the rock acts to create more surface energy in the medium (Bercovici & Ricard 2005). This theory avoids the rather nebulous approach of assuming that damage exists as some state variable without discerning how the damage is manifested in the system. An important question to ask though is how these manifestations of damage are related to the observed microphysics of damage. Void generation or Mode I cracks are well observed in brittle behaviour (and possibly brittle-ductile behaviour) at low pressures (Kohlstedt *et al.* 1995). At greater depths (~15–18 km) the overburden pressure acts to suppress rock dilation (Kohlstedt *et al.* 1995), so the relevance of void generating damage to mid to deep lithospheric depths is of questionable importance (unless these voids were somehow filled with fluid which seems difficult given the petrologically inferred lack of water in the lithosphere). The relevance of grain size (and specifically the reduction in grain size) to the rheological properties of rocks and shear localization is well demonstrated (Karato 1983; Jin *et al.* 1998). The cause of reduction in grain size is generally associated with dynamic recrystallization (Karato *et al.* 1980). The relation between deformational work partitioning into grain size reduction in the two-phase damage theory and the microphysics of dynamic recrystallization is not currently understood and remains an important and essential consideration for understanding this type of damage (Bercovici & Ricard 2005). Two-phase damage theory was recently incorporated by Bercovici & Ricard (2005) into a simple source-sink flow to determine the efficacy of both void and fineness generating damage to produce plate like flow (solid body translation, toroidal motion) in a 2-D sheet. Their results showed that fineness generating damage was very successful at generating plate like flow, while void generating damage merely enhanced the dipolar source-sink flow field and was, therefore, highly unplate-like (Bercovici & Ricard 2005).

In this study, we will extend the previous work done on two-phase damage theory (Bercovici *et al.* 2001b; Ricard *et al.* 2001; Ricard & Bercovici 2003; Bercovici & Ricard 2003, 2005) by incorporating two-phase damage physics into a coupled plate-mantle model that is convectively driven. Our model formulation attempts

to add another level of sophistication to previous work, but is still obviously a simplification of fully variable viscosity convection and hence convection in Earth. Our goal though is to ascertain the general plate behaviours (i.e. narrow boundary zones separated by rigid plate interiors) of two-phase damage theory in a convective system were we assume that the rheology is controlled by a single strong layer in the lithosphere. A simple but important addition to the evolution equation for fineness (inverse grain size) is added to the theory by the incorporation of a healing term. Healing in the fineness equation is representative of coarsening or grain growth in a poly-crystalline rock and is a well studied phenomena in Earth materials and in general (Karato 1989). One expects that increasing the fraction of deformational work which goes into generating damage as well as increasing the sensitivity of viscosity to grain size will allow for greater shear localization. In contrast, we expect that increasing the healing parameter will result in moderating the amount of localization within the lithosphere but potentially allow for a strengthening of the plate interior with respect to the boundaries. While we do find that increasing the sensitivity of viscosity to grain size results in greater shear localization, the effect of increasing the healing parameter and fraction of deformational work stored as damage exhibits some behaviours different than expected. We also develop a series of scaling analyses for how the convergent and divergent boundary widths depend upon the model parameters and what damage inputs control the boundary behaviour. Scaling analyses are useful guides for exploring parameter space and understanding the different localization behaviours demonstrated by the numerical experiments. In this study, we will inspect the ability of both void and fineness generating damage in a convectively driven system to produce Earth-like features and determine what implications this may have for plate tectonics on Earth.

## 2 MODEL FORMULATION

The goal of this paper is to determine the effectiveness of different damage formulations in producing plate-like flow in a 2-D convectively driven system; this will involve determining how convergent and divergent zones as well as slowly deforming ‘plate-like’ zones respond to the given damage parameters. To that end we will employ a lithosphere-mantle coupling model whereby a Newtonian mantle is overlain by a two-phase damage rheology lithospheric layer (Fig. 1); this model is similar to the formulation of Weinstein & Olson (1992). The lithospheric layer is of constant thickness *h* and is much smaller than the thickness of the underlying mantle *D*, ($h \ll D$). Material which flows across the boundary between

<figure><img src="image_2.png" /><figcaption>**Figure 1.** The thin sheet model employed in this study along with the vertical viscosity profile.</figcaption></figure>
the lithosphere and mantle undergoes an instantaneous rheological transformation to the rheological properties of the region into which it flows. The boundary between the two regions is a no-slip surface such that stresses transmitted from the Newtonian layer to the lithosphere drive the flow in the surface layer. In the context of the thin sheet approximation that will be described hereafter, the temperature at the surface of the Newtonian mantle is assumed to be uniform and the thin layer is isothermal. With the assumed thermal profile of the lithosphere it follows that the top layer has no buoyancy variations to drive flow in the model. Contrary to similar model formulations (Weinstein & Olson 1992; Weinstein 1996) the lithosphere is recycled back into the mantle by calculating the vertical velocity at the plate-mantle interface, this avoids the rather unphysical situation where lithospheric material would build up at a subduction zone and become depleted at a ridge.

## 2.1 Mantle convection model

In order to focus on the effect of the various forms of damage on convection we employ a rather simple model for mantle convection, that is isoviscous, Boussinesq, infinite Prandtl number, Rayleigh-Benard convection. The equations for mass, momentum and energy are

$$ \nabla \cdot \underline{\mathbf{v}}_m = 0 \qquad (1) $$

$$ 0 = -\nabla P_m + \nabla^2 \underline{\mathbf{v}}_m + Ra\theta\hat{\mathbf{z}} \qquad (2) $$

$$ \frac{\partial \theta}{\partial t} + \underline{\mathbf{v}}_m \cdot \nabla \theta - w_m = \nabla^2 \theta, \qquad (3) $$

where the subscript *m* denotes variables in the mantle, and $P_m$, $\underline{\mathbf{v}}_m$, $w_m$ and $\theta$ are the mantle non-hydrostatic pressure, velocity, vertical velocity and temperature perturbation, respectively (Schubert et al. 2001). The constant $Ra$ is the Rayleigh number for the system. The convection equations above have also been non-dimensionalized by $x = Dx'$, $t = (D^2/\kappa)t'$, $\underline{\mathbf{v}}_m = (\kappa/D)\underline{\mathbf{v}}_m'$ and $P_m = (\mu_m \kappa/D^2)P_m'$, and the primes have subsequently been dropped. The constants $D$, $\kappa$ and $\mu_m$ are the depth of the convecting system, thermal diffusivity and mantle viscosity, respectively. The above equations with the applied boundary conditions (to be discussed later) determine the behaviour of the Newtonian mantle in our model.

## 2.2 Two-phase damage theory: review and current formulation

The two-phase damage equations originate from a series of papers (Bercovici et al. 2001a,b; Ricard et al. 2001), with subsequent papers refining various aspects of the theory (Bercovici & Ricard 2003, 2005; Ricard & Bercovici 2003). The equations presented below are in the geologically applicable 'void limit' as discussed in Ricard & Bercovici (2003), whereby the void phase has zero density, pressure and viscosity.

### 2.2.1 Mass

The equation for mass conservation is

$$ \frac{D\phi}{Dt} = \frac{\partial\phi}{\partial t} + \mathbf{v} \cdot \nabla\phi = (1-\phi)\nabla \cdot \mathbf{v}, \qquad (4) $$

where $\mathbf{v}$ is the matrix (or rock phase) velocity and $\phi$ is the porosity. In the void limit the velocity of the void (or secondary) phase is

the same as the matrix velocity, hence we do not need to solve independent mass and momentum equations for the void phase.

### 2.2.2 Momentum

The momentum equations for the matrix are

$$ 0 = -\nabla[(1-\phi)P - \sigma\alpha] + \nabla \cdot [(1-\phi)\underline{\mathbf{\tau}}] - (1-\phi)\rho g \hat{\mathbf{z}}, \qquad (5) $$

where $P$ is the matrix pressure, $\sigma$ is the surface tension, $\alpha$ is the interfacial area density, $\underline{\mathbf{\tau}}$ is the deviatoric matrix stress and $\rho$ is the density of the matrix phase. The deviatoric stress is given by

$$ \underline{\mathbf{\tau}} = \mu_L \left[ \nabla \mathbf{v} + \left[ \nabla \mathbf{v} \right]' - \frac{2}{3} \left( \nabla \cdot \mathbf{v} \right) \underline{\mathbf{I}} \right], \qquad (6) $$

where $\underline{\mathbf{I}}$ is the identity matrix and $\mu_L$ is the top layer viscosity. The interfacial area density has been suggested to go as

$$ \alpha = A\eta(\phi), \quad \eta(\phi) = \phi^a(1-\phi)^b, \qquad (7) $$

where $A$ is the inverse void or grain size (depending on the curvature of the interface), $\eta(\phi)$ is a dimensionless function of porosity and $a, b$ are constants $\le 1$. The above equation satisfies the requirement that interfacial area goes to zero at the different limits of porosity, and simple micromechanical models have shown how $A$ relates to inverse grain/void size (Bercovici et al. 2001a). Following the previous formulations of the matrix rheology we assume the lithospheric viscosity is given by

$$ \mu_L = \mu_o \left( \frac{A_{\text{ref}}}{A} \right)^m, \qquad (8) $$

where $A_{\text{ref}}$ is a reference value for fineness and $\mu_o$ is reference viscosity of the lithosphere. The viscosity exponent $m$ is a dimensionless positive constant, and here we consider a reasonable range $m$ (i.e. $1 \le m \le 3$) assuming a grain size sensitive deformation mechanism (e.g. diffusion creep) is in operation. Given the above momentum equations (eq. 5) we find that the effective matrix viscosity is given by

$$ \mu_{\text{eff}} = (1-\phi)\mu_L = \mu_o(1-\phi)\left(\frac{A_{\text{ref}}}{A}\right)^m, \qquad (9) $$

where the $(1-\phi)$ factor in the viscosity essentially arises from the deviatoric stress term in (eq. 5).

### 2.2.3 Energy

The energy equation is separated into two coupled equations which govern the evolution of thermal energy and the rate of work done on the interface by pressure, surface tension and deformational work (Bercovici & Ricard 2005). The evolution of thermal energy is

$$ (1-\phi)\rho c \frac{\mathrm{D}T}{\mathrm{D}t} - T \frac{\mathrm{D}}{\mathrm{D}t} \left( \alpha \frac{\mathrm{d}\sigma}{\mathrm{d}T} \right) - T\alpha \frac{\mathrm{d}\sigma}{\mathrm{d}T} \nabla \cdot \mathbf{v} = \\ Q - \nabla \cdot \mathbf{q} + B \left( \frac{\mathrm{D}\phi}{\mathrm{D}t} \right)^2 + \sigma\eta k_A A^p + (1-f)\Psi, \qquad (10) $$

where $T$ is the temperature, $c$ is the matrix heat capacity, $-b\sigma/\mathrm{D}T$ is the interfacial entropy per unit area, $\mathbf{q}$ is an energy flux vector (e.g. heat diffusion) and $Q$ is an intrinsic heat source. The quantity $B$ is positive, has units of viscosity, and is related to the bulk viscosity term in the two-phase theory of McKenzie (1984); the term proportional to $B$ represents irreversible work done during
isotropic compression or dilation. Following previous formulations for B we arrive at

$$B = \frac{K \mu_o}{\phi(1 - \phi)}, \qquad (11)$$

where K is a dimensionless factor accounting for pore or grain geometry and is typically $O(1)$ (Bercovici et al. 2001a; Bercovici & Ricard 2005). The $(1 - \phi)$ in denominator of (eq. 11) may not be necessary in the void limit, but the range of porosities in our study make this term negligible anyways. The quantity $k_A$ represents the rate of grain growth (Karato 1989) and $\eta(\phi)$ is a dimensionless function of porosity (see eq. 7); the term proportional to $k_A$ represents the contribution of irreversible loss of interfacial area (via grain growth) to the evolution of thermal energy. The exponent p in the term proportional to $k_A$ is related to surface tension driven grain growth which will be further elaborated on later in this section. The viscous deformational work is given by

$$\Psi = (1 - \phi)\nabla v : \tau, \qquad (12)$$

a fraction $1 - f$ of which is partitioned into dissipative heating. The evolution of energy associated with work done on the interface and interface generation is

$$\sigma \frac{D\alpha}{Dt} = -P \frac{D\phi}{Dt} - B \left(\frac{D\phi}{Dt}\right)^2 - \sigma\eta k_A A^p + f\Psi. \qquad (13)$$

Since changes in $\alpha$ can be manifested as either changes in porosity or changes in inverse grain/void size we can decompose (eq. 13) into one damage evolution equation for porosity and one for grain size (Bercovici & Ricard 2005). The damage equation for porosity is

$$\sigma A \frac{d\eta}{d\phi} \left(\frac{D\phi}{Dt}\right) = -P \left(\frac{D\phi}{Dt}\right) - B \left(\frac{D\phi}{Dt}\right)^2 + f_{\phi} \Psi \qquad (14)$$

and the evolution equation for fineness is

$$\frac{DA}{Dt} = \frac{f_A}{\sigma\eta}\Psi - k_A A^p, \qquad (15)$$

where $f = f_{\phi} + f_A$ must be less than one. The association of the term proportional to $k_A$ with surface tension driven grain growth (and hence the curvature of grains) implies a direct relationship to the size of the grains, hence we assume the $k_A$ healing term is only in the equation governing the evolution of grain size (eq. 15) and not in eq. (14). The term proportional to $k_A$ in (eq. 15) acts to decrease the fineness (or increase grain size), and decreasing fineness leads to an increase in viscosity and, therefore, heal the damaged zone (i.e. plate boundary). We will consider cases where both the void- and fineness-generating damage mechanisms operate by them self as well as in tandem with each other. The partitioning fraction of damage in the void-generating case is given by

$$f_{\phi} = f^* \frac{(D\phi/Dt)^2}{\gamma + (D\phi/Dt)^2}, \qquad (16)$$

where $f^*$ is the maximum permissible $f_{\phi}$, $\gamma$ controls the variability of $f_{\phi}$ and $f_{\phi}$ is assumed to depend on an even power of $D\phi/Dt$ since it must be positive definite. As previously discussed in Bercovici et al. (2001b) and Bercovici & Ricard (2003), the above equation for $f_{\phi}$ precludes singular solutions of the porosity damage equation (eq. 14) in areas of zero void growth (i.e. $D\phi/Dt = 0$).

The fineness generating damage equation presented in this paper is different from its previous form (Bercovici & Ricard 2005) in that it now explicitly contains a healing term. In the case where fineness represents inverse grain size (the only case considered in

this paper) the healing term would be representative of coarsening or grain growth (Karato 1989). Whether A represents inverse grain or void size depends upon the values of a and b in (eq. 7) as discussed previously. Similar to Bercovici & Ricard (2005) we find that variations in these parameters do not significantly affect the outcome of our results, and making the values of a smaller (which puts the fineness parameter in the inverse grain size regime) are easily offset by increasing $f_A$. We will, therefore, associate the fineness parameter with inverse grain size in this paper and attempt to use constraints from grain size evolution to determine parameters in our evolution equation.

Previous analyses that have incorporated a damage evolution equation similar to the fineness evolution equation suggested the association of this equation to inverse grain size as well (Tackley 2000b; Auth et al. 2003), and in both cases choose $p = 1$ in eq. (15). The fineness evolution equation (eq. 15) in the absence of damage satisfies

$$\frac{DA}{Dt} \approx -k_A A^p \Rightarrow \frac{1}{A^{p-1}} \approx (p-1)k_A t. \qquad (17)$$

Since we want to associate A with inverse grain size, $A = 1/g$, and the evolution of grain size undergoing grain growth is often found to vary as $(k_A t)^{\frac{1}{p-1}}$ (e.g. Karato 1989; Evans et al. 2001) we therefore, choose $p = 3$. Choosing $p = 1$ leads to a logarithmic integration where the grain growth rate increases with time, while choosing $p = 3, 4$ allows for grain growth to be driven by surface tension and curvature. This simple analysis will provide the guidance for our choice of the healing exponent in eq. (15). Therefore, the pre-factor $k_A$ plays the role of the growth rate, which is experimentally found to be a function of temperature (Karato 1989). As we will discuss in the next section, temperature variations in our model lithosphere are not included which implies that the growth rate $k_A$ in our model is constant. The lithosphere is likely strongest at intermediate depths (30–70 km at 100 Myr) and, therefore, choosing a grain growth healing rate associated with this depth would be most appropriate in the thin-sheet formulation. The first term in eq. (15) allows for the input of deformational work to go into the reduction of grain size and increase the surface energy in the system. Grain size reduction in Earth minerals is experimentally found to occur by dynamic recrystallization (Karato et al. 1980), but there exists no complete development for the relationship between grain size and these processes (Hall & Parmentier 2003). Our evolution equation takes a similar form as previous studies, but we caution that our grain size reduction term comes about from the partitioning of deformational work between reversible and irreversible processes and not from a model that is meant to specifically represent any one experimentally determined grain size reduction mechanism. Recent work by Austin & Evans (2007) has also suggested that grain size reduction is determined by the rate of deformational work, and their proposed evolution equation for grain size (which is very similar to eq. 15) reproduces experimental results quite well.

## 2.3 Lithosphere model formulation

The thin layer that represents the lithosphere is assumed to have a free-slip, impermeable surface on top and a no-slip interface with the Newtonian layer. Since our lithosphere thickness is assumed small we will employ the thin-sheet approximation to solve for the dynamics of this layer (England & McKenzie 1982; Wdowinski et al. 1989; Ribe 1992; Weinstein & Olson 1992; Lemery et al. 2000). In the thin-sheet approximation, vertical gradients in horizontal velocity are assumed negligible in comparison to the horizontal
gradients. The thin-sheet formulation, therefore, involves determining the lithosphere's vertically averaged horizontal velocity and the vertical velocity at the interface with the Newtonian layer. We thus examine the horizontal and vertical components of the momentum equation in a vertically averaged sense. The x and z components of the momentum equation in the lithosphere are:

$$0 = \frac{\partial}{\partial x} [\sigma\alpha + (1-\phi)(\tau_{xxL} - P)] + \frac{\partial}{\partial z} [(1-\phi)\tau_{xzL}] \quad (18)$$

$$0 = \frac{\partial}{\partial z} [\sigma\alpha + (1-\phi)(\tau_{zzL} - P)] + \frac{\partial}{\partial x} [(1-\phi)\tau_{xzL}], \quad (19)$$

where the subscript L designates that the stresses are in the thin layer and we have neglected the buoyancy force of the top layer. In the $\hat{z}$ eq. (19) we will neglect the shear stress term since it goes as $O(h/L)$ (since $\partial\tau_{xz}/\partial x \sim h u_L/L$, where $u_L$ is the horizontal velocity in the lithosphere), which implies $\sigma_{zzL}$ is constant. Taking the vertical average of both equations (and assuming the lithospheric thickness $h$ is constant) results in

$$0 = \frac{\partial}{\partial x} [\sigma \bar{\alpha} + (1 - \bar{\phi})(\bar{\tau}_{xxL} - \bar{P})] - \frac{1}{h} [(1 - \bar{\phi})\tau_{xzm}|_{z=D}] \quad (20)$$

$$\begin{aligned} \frac{\partial \sigma_{zzL}}{\partial z} &= 0 \Rightarrow \sigma_{zzL}|_{D}^{D+h} = 0 \Rightarrow [\sigma \bar{\alpha} + (1 - \bar{\phi})(\bar{\tau}_{zzL} - \bar{P})] \\ &= \tau_{zzm}|_{z=D} \end{aligned} \quad (21)$$

using the continuity of vertical and shear stress at the boundary between the lithosphere and Newtonian layer. We define $\sigma_{zzL} = \sigma\bar{\alpha} + (1-\bar{\phi})(\bar{\tau}_{zzL}-\bar{P})$ to be the total vertical stress in the lithosphere, and the subscript $m$ refers to stresses in the Newtonian mantle. Eq. (21) results from our approximation that $\sigma_{zzL}$ is constant in the lithosphere, hence the vertically averaged $\sigma_{zzL}$ is equal to $\sigma_{zzL}$ at $z = D+h$ (and also at $z = D$). To close the system of equations we also make the approximation that the average of a product is the product of the averages (e.g. $\phi\bar{\tau}_L = \bar{\phi}\bar{\tau}_L$). The vertically averaged stress is given by integrating eq. (6), and making use of the assumption that the vertical velocity is zero at the top boundary; this results in

$$\bar{\tau}_{xxL} = \bar{\mu}_L \left[ \frac{4}{3} \frac{\partial \bar{u}_L}{\partial x} + \frac{2}{3} \frac{w_L}{h} \Big|_{z=D} \right] \quad (22)$$

$$\bar{\tau}_{zzL} = \bar{\mu}_L \left[ -\frac{4}{3} \frac{w_L}{h} \Big|_{z=D} - \frac{2}{3} \frac{\partial \bar{u}_L}{\partial x} \right]. \quad (23)$$

The vertically averaged evolution equations for damage are

$$\bar{P} = -\sigma \bar{A} \frac{d\bar{\eta}}{d\bar{\phi}} - \frac{K \mu_o}{\bar{\phi}} \nabla \cdot \bar{\mathbf{v}} + F, \quad (24)$$

where

$$\nabla \cdot \bar{\mathbf{v}} = \nabla_H \cdot \bar{\mathbf{v}} - \frac{w_L}{h} \Big|_{z=D}, \quad \bar{\mathbf{v}} = (u_L, 0, w_L) \quad (25)$$

and

$$F = f^* \frac{(1 - \bar{\phi})^2 \nabla \cdot \bar{\mathbf{v}}}{\gamma + (1 - \bar{\phi})^2 (\nabla \cdot \bar{\mathbf{v}})^2} \bar{\Psi} \quad (26)$$

under the same assumptions that lead to eqs (20) and (21). Substituting in the expressions for $\bar{\tau}_{xxL}$, $\bar{\tau}_{zzL}$ and $\bar{P}$ eqs (22–24) into the momentum eqs (20 and 21) we get the following system of equations for $\frac{\partial u_L}{\partial x}$ and $\frac{w_L}{h}|_{z=D}$:

$$\begin{pmatrix} a_m & b_m \\ -b_m & -a_m \end{pmatrix} \begin{pmatrix} \frac{\partial u_L}{\partial x} \\ \frac{w_L}{h} \end{pmatrix} = \begin{pmatrix} \frac{1}{h} \int_0^x (1-\bar{\phi}) \tau_{xzm} |_{z=D} dx' + \tau_{xz_o} - \sigma \bar{A} \lambda \\ \tau_{zzm} |_{z=D} - \sigma \bar{A} \lambda \end{pmatrix} + \begin{pmatrix} F \\ F \end{pmatrix}, \quad (27)$$

where the constant $\tau_{xz_o}$ is determined using periodic boundary conditions, and

$$\begin{aligned} a_m &= (1 - \bar{\phi}) \left[ \frac{4}{3} \mu_L + \frac{K \mu_o}{\bar{\phi}} \right] \\ b_m &= (1 - \bar{\phi}) \left[ \frac{2}{3} \mu_L - \frac{K \mu_o}{\bar{\phi}} \right] \\ \lambda &= \left[ \bar{\eta} + (1 - \bar{\phi}) \frac{d\bar{\eta}}{d\bar{\phi}} \right]. \end{aligned} \quad (28)$$

We next consider the non-dimensionalization of the equations governing the behaviour of the lithosphere. We will employ the same dimensional scales as used for the mantle, except that we will include an additional length scale $A_{ref}$ (m$^{-1}$) to non-dimensionalize the fineness (inverse grain size) in order to keep the macroscopic length scale (D) separate from the microstructural one since they differ by about nine orders of magnitude. Substituting the aforementioned scales into our governing equations (with the bar designating a vertically averaged quantity subsequently dropped) results in the following non-dimensional equations

$$\begin{pmatrix} a_m & b_m \\ -b_m & -a_m \end{pmatrix} \begin{pmatrix} \frac{\partial u_L}{\partial x} \\ \frac{w_L}{h} \end{pmatrix}|_{z=1} = \begin{pmatrix} \frac{1}{h} \int_0^x (1-\phi) \tau_{xzm} |_{z=D} dx' + \tau_{xz_o} - \hat{\sigma} A \lambda \\ \tau_{zzm} |_{z=1} - \hat{\sigma} A \lambda \end{pmatrix} + \mu_R \begin{pmatrix} F \\ F \end{pmatrix}, \quad (29)$$

where

$$\begin{aligned} a_m &= (1 - \phi) \left[ \frac{4}{3} \frac{\mu_R}{A^m} + \frac{K \mu_R}{\phi} \right], \\ b_m &= (1 - \phi) \left[ \frac{2}{3} \frac{\mu_R}{A^m} - \frac{K \mu_R}{\phi} \right]. \end{aligned} \quad (30)$$

We have defined the following non-dimensional quantities in the above equations as

$$\hat{\sigma} = \frac{\sigma A_{\text{ref}} D^2}{\mu_m \kappa}, \quad (31)$$

$$\mu_R = \frac{\mu_o}{\mu_m}, \quad (32)$$

where $\mu_o$ is the reference lithospheric viscosity. The non-dimensional version of the mass equation is exactly the same as the dimensional case (eq. 4), while the fineness evolution equation becomes

$$\frac{DA}{Dt} = \frac{f_A \mu_R}{\eta \hat{\sigma}} A^{-m} \Psi^* - \hat{k}_A A^p, \quad (33)$$

where

$$\hat{k}_A = \frac{k_A A_{\text{ref}}^{p-1} D^2}{\kappa}, \quad (34)$$

$$\Psi^* = (1 - \phi) \left\{ \nabla \mathbf{v} : [\nabla \mathbf{v} + [\nabla \mathbf{v}]'] - \frac{2}{3} (\nabla \cdot \mathbf{v})] \right\}. \quad (35)$$

The general scheme for solving the model system proceeds along the following steps. (1) Solve for the velocity and stress in the Newtonian mantle (eqs 1–2) given the current density distribution subject to free-slip bottom and no-slip top boundary conditions. (2) The flow calculated in (1) generates a stress at the lithosphere–mantle interface that is used to calculate the velocity in the lithosphere via (eq. 29). (3) The lithospheric velocity is then used to drive cavity flow in the Newtonian mantle, and the velocity and stress in the Newtonian mantle are the sum of solutions from (1) and (3). (4) We
iterate on (2) and (3) until the L2 norm of the difference between successive iterations reaches $10^{-6}$ or less, at which point the lithosphere and mantle are assumed to be in mechanical equilibrium. (E) Finally we update the temperature field in the mantle (eq. 3) and the damage variables in the lithosphere (eqs 4 and 33) with the equilibrium velocities determined at the end of (4). Convection in the Newtonian mantle is solved numerically via a spectral method that employs a Propagator Matrix method for solving Stokes flow (Hager & O'Connell 1981) and finite difference for the temperature equation. The cavity flow velocity field is calculated using a Propagator Matrix method as well. The sensitivity of lithospheric and convective dynamics to variations in damage parameters we will consider in this study include varying $f_\phi$ and $f_A$, $k_A$ and $m$, which are the fraction of deformational work partitioned into generating voids and fineness, grain growth healing rate, and the viscosity exponent for grain size sensitivity, respectively.

## 2.4 Scaling analysis for boundary width/strain-rate

In the study of plate generation, an important goal is the formation of narrow zones of weakness (plate boundaries) which separate strong plate-like interiors, especially when the analysis is confined to two dimensions. Before proceeding with a large number of numerical experiments it is useful to develop a scaling analysis that can guide our exploration of parameter space as well as provide insight into the different behaviours our numerical experiments produce. We are interested in determining how the plate boundary widths depend upon the free parameters in our model at steady state; to that end, we look at the fineness evolution eq. (33) in the absence of time-dependence:

$$u_L \frac{dA}{dx} = \frac{f_A \mu_R}{\hat{\sigma} \eta} A^{-m} \dot{\epsilon}^2 - \hat{k}_A A^p. \qquad (36)$$

We assume functions with the proper symmetries for the velocity and fineness in the vicinity of the plate boundary, and then substitute them into eq. (15) to solve for the boundary widths. We thus assume

$$\begin{aligned} u_L &= \pm u_{\max} \tanh(x/\delta) \\ A &= A_o \exp(-x^2/\delta^2), \end{aligned} \qquad (37)$$

where the $\pm$ in the velocity equation distinguishes the relations between convergent (-) and divergent (+) zones, $u_L$ is the surface velocity and $A_o$ is the steady-state maximum amplitude of fineness. The constant $\delta$ is the characteristic width for each boundary (convergent and divergent) and is assumed to be equal for both the fineness and velocity field (i.e. $\delta_A = \delta_v$). Multiplying eq. (36) through by $A^m$ and substituting the above functional forms (eq. 37) into eq. (36) gives

$$\begin{aligned} & \pm \frac{2u_{\max}}{\delta^2} A_o^{m+1} x \tanh(x/\delta) \exp[-(m+1)x^2/\delta^2] \\ & = \frac{f_A \mu_R u_{\max}^2}{\hat{\sigma} \eta \delta^2} \cosh^{-4}(x/\delta) - \hat{k}_A A_o^{m+p} \exp[-(m+p)x^2/\delta^2]. \end{aligned} \qquad (38)$$

We next integrate eq. (38) around the plate boundary over a domain much greater than the boundary width $\delta$. We make the assumption that the integral of eq. (38) from $[-\delta/2 - L, +\delta/2 + L]$ is approximately equal to the integral with the limits $(-\infty, +\infty)$ to facilitate evaluation of the above integral (i.e. $\delta \ll L$). We also approximate $\tanh(x/\delta)$ by $x/\delta$ which doesn't produce any significant error since the Gaussian in the relevant term decays much faster than the linear term grows. The resultant integration of eq. (38) and organization

of terms yields

$$\begin{aligned} & \sqrt{\frac{\pi}{m+p}} \hat{k}_A A_o^{m+p} + (\mp) \frac{\sqrt{\pi}}{(m+1)^{3/2}} A_o^{m+1} \left(\frac{u_{\max}}{\delta}\right) \\ & - \frac{4}{3} \mu_R \frac{f_A}{\hat{\sigma}\eta} \left(\frac{u_{\max}}{\delta}\right)^2 = 0, \end{aligned} \qquad (39)$$

with which we can estimate $\delta$ (or $u_{\max}/\delta$, the boundary strain-rate) in terms of our governing free parameters. The term proportional to $(u_{\max}/\delta)^0$ in eq. (39) is associated with healing, the term proportional to $(u_{\max}/\delta)^1$ is associated with advection $(-)$ is for divergent and $(+)$ is for convergent], and the term proportional to $(u_{\max}/\delta)^2$ is associated with deformational work. Solving for the boundary widths for convergence ($\delta_c$) and divergence ($\delta_d$) gives

$$\delta_{d,c} = \frac{\pm 1 + \sqrt{1 + \frac{16}{3\sqrt{\pi}} \frac{(m+1)^3}{(m+p)^{1/2}} A_o^{p-m-2} \mu_R \frac{f_A}{\hat{\sigma}\eta} \hat{k}_A}}{2 \frac{(m+1)^{3/2}}{(m+p)^{1/2}} A_o^{p-1} \frac{\hat{k}_A}{u_{\max}}}. \qquad (40)$$

It is clear that eq. (40) provides a scaling law for both boundary width ($\delta$) and boundary strain-rate $u_{\max}/\delta$. We focus our analysis on both boundary widths and strain-rates for each numerical experiment to characterize the degree and type of localization present. The scaling analysis for $\delta$ above depends upon the steady-state fineness amplitude ($A_o$), and $A_o$ is clearly a function of the input damage parameters.

## 2.5 Scaling analysis for fineness magnitude

Without including a healing term in the damage eq. (15) it is obvious that fineness would continue to increase (i.e. grain size continues to decrease) without bound which is an unphysical scenario. Given that grain growth acts to increase the average grain size, even in the presence of deformation (Karato 1983, 1989) it is important to consider its effects on steady state fineness. We first examine how the magnitude of fineness depends upon the damage parameters ($f_A$, $\hat{k}_A$ and $m$), as well as consider the stability of the fineness evolution eq. (15) to an infinitesimal perturbation. We seek the maximum value of fineness at the boundaries where $dA/dx$ is zero, and advection does not affect the fineness magnitude. Expanding the fineness field around its steady state value ($A_o$)

$$A = A_o + \epsilon A'(t), \qquad (41)$$

where $A'$ is the perturbation fineness and $\epsilon \ll 1$, and substituting this into

$$\frac{dA}{dt} = f_A A^m \frac{\tau^2}{\hat{\sigma} \eta \mu_R} - \hat{k}_A A^p \qquad (42)$$

we determine how the stability of steady state fineness is governed by the damage parameters. To 0th order in $\epsilon$ [assuming that the stress ($\tau$) is constant] the steady-state fineness goes as

$$A_o = \left( \frac{\tau^2}{\hat{\sigma} \eta \mu_R} \right)^{\frac{1}{p-m}} \left( \frac{f_A}{\hat{k}_A} \right)^{\frac{1}{p-m}}, \qquad (43)$$

and to first order in $\epsilon$

$$\frac{dA'}{dt} = (m-p)\hat{k}_A A_o^{p-1} A'. \qquad (44)$$

Two interesting points can be made. First, the steady-state fineness is predicted to be a function of the ratio of $f_A$ and $\hat{k}_A$, second, stability of the steady-state to perturbations is determined solely by the sign of $m-p$. The scaling above for the magnitude of the steady-state fineness also provides a diagnostic to determine whether or not fineness magnitudes are determined by the competition between deformational work and healing.
## 2.6 A note on boundary widths and the influence of temperature

Some important points can be made regarding the results of the scaling analysis for plate boundary widths. Convergent boundaries are predicted to be smaller than divergent boundaries for a given plate velocity, and this difference arises predominately from the difference in sign of the velocity field (see eq. 37). Consideration of data from global seismicity would imply the opposite of the above prediction, specifically that divergent boundaries are narrower than convergent boundaries (Dumoulin et al. 1998). The difference in boundary widths between divergent and convergent zones seen on Earth are likely to not arise solely due to the difference in sign of the divergence field, but it also points out the importance of temperature variations in the lithosphere. The characteristic temperature of convergent and divergent boundaries in the lithosphere is very different, and this necessarily plays an important role for the variations in boundary widths between convergence and divergence. The effects of horizontal variations in temperature within the plate and their subsequent influence on viscosity are not explicitly included in this model due to the assumed isothermal structure of the plate. While we include the zeroth order influence of temperature dependent viscosity in that $\mu_o > \mu_m$ the existence of horizontal variations in viscosity due to horizontal variations in temperature likely influences the difference between our scaling law predictions and Earth. Inspection of the scaling law in eq. (40) suggests that allowing for temperature variations in the lithosphere would most readily affect the viscosity ($\mu_o$) and the healing rate ($k_A$); an increase (decrease) in temperature would decrease (increase) viscosity while increasing (decreasing) healing rate. The divergent boundary width with temperature variations would, therefore, decrease compared to the isothermal case and vice versa for convergent boundaries. Given that the scaling analysis is derived by considering local variations we could just as easily assume that $\mu_o$ and $k_A$ are different for convergent and divergent boundaries. The scaling law, therefore, seems capable of predicting boundary width variations seen on Earth when the appropriate lithospheric temperature variations are taken into consideration. We will discuss more later about the influence of variations in viscosity due to temperature and its inclusion in subsequent models.

## 3 NUMERICAL RESULTS

We next explore numerical solutions of our model in a sampling of parameter space. All of our convection simulations have a simple Rayleigh-Bernard heating mode, and the Rayleigh number is $10^6$ with an aspect ratio of two for each case. We choose to use this simple bottom heating mode in order to take advantage of the symmetries and steady-state behaviour of this convective system. In all cases we initiated our damage rheology calculations with a Newtonian convection simulation that has reached steady-state. For each set of parameters investigated (see Table 1 for list of all parameters), we attempted to find steady-state solutions to find robust measurements of the plate characteristics (boundary widths, concentration of deformation and rigidity of plate interiors) that we are most interested in understanding. We also choose the lithospheric plate viscosity to be four orders of magnitude greater than the viscosity of the Newtonian mantle (i.e. $\mu_R = 10^4$). With typical values of mantle scales ($\mu_m = 10^{21}$ Pa s, $\sigma = 1$ N m$^{-1}$, $\kappa = 10^{-6}$ m$^2$ s$^{-1}$, $D \sim 10^6$ m and $A_{\text{ref}} \sim 10^3$ m$^{-1}$ implying a grain size of 1 mm) the non-dimensional surface tension, $\hat{\sigma}$ is approximately 1, and we will assume that it is throughout the calculations. When

<table><caption>**Table 1.** List of the values for the model parameters: $f_\phi$ is the fraction of deformational work that goes into void generating damage, $f_A$ is the fraction of deformational work that goes into fineness generating damage, $m$ controls the sensitivity of viscosity to fineness, $\hat{k}_A$ is the healing rate related to grain growth and $\gamma$ modulates the amount of void generating damage.</caption><thead><tr><th>f<sub>φ</sub></th><th>f<sub>A</sub>(×10<sup>-4</sup>)</th><th>m</th><th>k&#x0302;<sub>A</sub></th><th>γ</th></tr></thead><tbody><tr><td>0</td><td>0.05</td><td>1.2</td><td>1</td><td>0</td></tr><tr><td>0</td><td>0.25</td><td>1.2</td><td>0.75, 1, 2.5</td><td>0</td></tr><tr><td>0</td><td>0.5</td><td>1.2</td><td>2, 2.5, 5, 10, 25, 50</td><td>0</td></tr><tr><td>0</td><td>0.75</td><td>1.2</td><td>2.5, 3, 3.75, 5, 7.5, 15</td><td>0</td></tr><tr><td>0</td><td>1</td><td>1.2</td><td>4</td><td>0</td></tr><tr><td>0</td><td>0.05</td><td>1.4</td><td>5</td><td>0</td></tr><tr><td>0</td><td>0.25</td><td>1.4</td><td>5</td><td>0</td></tr><tr><td>0</td><td>0.5</td><td>1.4</td><td>10, 25, 50, 100</td><td>0</td></tr><tr><td>0</td><td>0.75</td><td>1.4</td><td>15</td><td>0</td></tr><tr><td>0</td><td>0.05</td><td>1.6</td><td>5</td><td>0</td></tr><tr><td>0</td><td>0.25</td><td>1.6</td><td>15, 37.5</td><td>0</td></tr><tr><td>0</td><td>0.5</td><td>1.6</td><td>25, 37.5, 50, 100</td><td>0</td></tr><tr><td>0</td><td>0.75</td><td>1.6</td><td>37.5</td><td>0</td></tr><tr><td>0</td><td>0.05</td><td>1.8</td><td>10</td><td>0</td></tr><tr><td>0</td><td>0.1</td><td>1.8</td><td>20</td><td>0</td></tr><tr><td>0</td><td>0.25</td><td>1.8</td><td>50</td><td>0</td></tr><tr><td>0</td><td>0.5</td><td>1.8</td><td>75, 80, 85, 100, 200</td><td>0</td></tr><tr><td>0</td><td>0.75</td><td>1.8</td><td>150</td><td>0</td></tr><tr><td>0</td><td>0.05</td><td>2</td><td>25</td><td>0</td></tr><tr><td>0</td><td>0.1</td><td>2</td><td>50</td><td>0</td></tr><tr><td>0</td><td>0.25</td><td>2</td><td>125</td><td>0</td></tr><tr><td>0</td><td>0.5</td><td>2</td><td>200, 250, 375</td><td>0</td></tr><tr><td>0</td><td>0.75</td><td>2</td><td>375</td><td>0</td></tr><tr><td>0</td><td>2.5</td><td>2</td><td>1250</td><td>0</td></tr><tr><td>0</td><td>25</td><td>2</td><td>12 500</td><td>0</td></tr><tr><td>0</td><td>50–1000</td><td>2</td><td>250 000</td><td>0</td></tr><tr><td>0</td><td>0.5</td><td>2.5</td><td>725, 750, 1000</td><td>0</td></tr><tr><td>0.5</td><td>0</td><td>0</td><td>0</td><td>10<sup>2</sup>, 1, 10<sup>-2</sup>, 10<sup>-4</sup></td></tr><tr><td>0.4</td><td>0.1</td><td>1.2, 2</td><td>2.5, 250</td><td>10<sup>-2</sup></td></tr><tr><td>0.25</td><td>0.25</td><td>1.2, 2</td><td>2.5, 250</td><td>10<sup>-2</sup>, 1, 10<sup>3</sup></td></tr><tr><td>0.2</td><td>0.3</td><td>1.2, 2</td><td>2.5, 250</td><td>10<sup>-2</sup></td></tr><tr><td>0.1</td><td>0.4</td><td>1.2, 2</td><td>2.5, 250</td><td>10<sup>-2</sup></td></tr></tbody></table>

the damage rheology is initiated the plate is assumed to have both constant porosity and fineness profiles, with values of 1 per cent and 1, respectively. The exponents (a and b) in the interfacial area density function $\eta(\phi)$ (7) are chosen to be 0.5 throughout the results section, and test cases with variations in a and b showed little to no difference in the system behaviour from the cases with $a = b = 0.5$. We employ 128 gridpoints in the horizontal direction and 50 gridpoints in the vertical direction. The chosen numerical resolution faithfully reproduced the results of test runs at higher resolutions but at significantly reduced computational times. Since it is difficult to glean information from numerical experiments that fail to converge to a stable solution our goal is to find a range of parameters that will produce numerically resolvable, steady-state behaviours, while still remaining in parameter regimes that produce Earth-like convective patterns. We will look at the following plate characteristics: divergent and convergent boundary width (we will measure this quantity by the full-width at half-max for the divergence field), as well as the concentration of deformation and rigidity of plate interiors as measured by the integrated diagnostic plateness [$Pl = (\pi - Pl^*)/(\pi - 2)$, where $Pl^* = (1/|u|_{\text{avg}}) \int |\dot{e}_{xx}| dx$] as defined by Weinstein & Olson (1992). The definition of plateness results in a sinusoidal velocity field having a plateness of 0, and boxcar velocity field having a plateness of 1. Looking at the individual
boundary widths will provide insight into how each deformational
zone responds to model parameters, while plateness will quantify
the narrowness of the boundaries as well as the rigidity of the plate
interiors.

**3.1 Void-generating damage**

The void-generating damage cases are distinguished by the variation
in the parameter $\gamma$ (and subsequently $f_{\phi}$), which primarily controls
the sensitivity of void-generating damage to the deformational work.
We examine cases for four different values of $\gamma$ (10$^2$, 1, 10$^{-2}$, 10$^{-4}$)
with constant $f^* = 0.5$, and show the results for $\gamma = 10^{-2}$ in Fig. 2.
Solutions for different $\gamma$ display similar trends: plateness is very low
($\sim$1, hence deformation is spatially distributed) and similar to the
case where there is no damage (i.e. $f_A = f_{\phi} = 0$), and porosity is
greater over ridges than in subduction zones (Fig. 2a). That voids are
generated over ridges and closed over subduction zones is indicative
of ridges and subduction zones undergoing tensional and compres-
sional environments, respectively. In the current void-generating
damage framework it seems that we are unable to generate narrow
plate boundaries and rigid plate interiors (Fig. 2b). Previous results
employing void-generating damage in simple-shear found that shear
localization occurred upon increasing $\gamma$ and deformational work in
tandem (Bercovici & Ricard 2003), but in the current study the de-
formational work is approximately constant since we don't vary the
Rayleigh number. Decreasing the magnitude of $\gamma$ results in a slight
increase in the rate at which voids are generated around divergent

<figure><img src="image_3.png" /></figure>

**Figure 2.** Numerical results for case with only void generating damage. In (a) we plot the porosity profile in the lithosphere, (b) we plot the surface velocity profile in the lithosphere and (c) we plot the temperature field. The parameter values are $f_{\phi} = 0.5$ and $\gamma = 10^{-2}$, and at $t = 0.05$ the results for surface velocity are in steady state. This example of void generating damage is representative of all such cases where the model results are similar to the case of a Newtonian rigid lid. The contour interval for the temperature field is 0.1 (and will be the same contour interval for all plots of the temperature field). The plateness (PI) for the velocity field is 0.1.

boundaries, but this variation produces no noticeable improvement
on localization for void-generating damage. This result corroborates
the work of (Bercovici & Ricard 2005) in that pure void-generating
damage is insufficient to allow for plate-like flow. It should be noted
that Bercovici & Ricard (2005) found that void-generating damage
was specifically unable to produce significant toroidal flow which is
not relevant to this study; however both their and our results suggest
that voids in the current two-phase framework do not allow for the
types of plate boundary behaviour seen on Earth. A possible way to
improve the plate-inducing capabilities of void-generating damage
would be to change the sensitivity of viscosity to porosity. This was
done in Bercovici & Ricard (2005) by specifying that the viscos-
ity undergoes a sharp drop when porosity goes above some critical
value, and Katz et al. (2006) employed a viscosity exponentially
dependent on porosity. While these different viscosity formulations
had varying degrees of success for the given model, we will not
explore variations on the porosity sensitivity to viscosity since it is
unclear of the relevance of the previous formulations to our study
and a more appropriate form is unknown at this time.

**3.2 Fineness-generating damage**

The fineness generating cases are distinguished by the variations
in the parameters *m*, *f*<sub>A</sub>, *k*<sub>A</sub> (see eq. 17). The healing exponent is
assumed constant at *p* = 3 throughout to reduce the number of
variables under consideration in this study. In the fineness evolution
eq. (33) the variables *f*<sub>A</sub> and *μ*<sub>R</sub> multiply each other and can be
treated as a single variable. Most of the results we present in the
following sections are for 0 ≤ *f*<sub>A</sub>*μ*<sub>R</sub> ≤ 1, and since *μ*<sub>R</sub> = 10<sup>4</sup>
this implies 0 ≤ *f*<sub>A</sub> ≤ 10<sup>-4</sup>. As we show in the following sections
the steady state results are not highly dependent upon the absolute
magnitude of *f*<sub>A</sub>*μ*<sub>R</sub>. Therefore, we can map out the general plate
generating behaviours of our model at smaller values of *f*<sub>A</sub>*μ*<sub>R</sub>, but
at reduced computational times compared to calculations at larger
*f*<sub>A</sub>*μ*<sub>R</sub>. Our results also suggest that *f*<sub>A</sub> is not required to be large
in order for significant amounts of localization to be generated.

3.2.1 *Variations in f*<sub>*A*</sub>

We consider variations in $f_A$ (which represents the partition of
work towards fineness inducing damage) at constant values of $\hat{k}_A$
(healing rate) and $m$ (viscosity exponent for grain size sensitiv-
ity). Increases in $f_A$ result in increases in the plate mobility (de-
fined as the ratio of plate velocity to rms velocity of the mantle)
and boundary strain-rate ($u_L/\delta_{d,c}$), while the boundary widths and
plateness remain fairly constant and insensitive to changes in $f_A$,
which is somewhat unexpected. These results imply that the fraction
of deformational work that goes into fineness-generating damage
is important for providing the reduction in strength (or increase
in strain-rate) necessary for plate boundary lubrication, but con-
tributes little to the localization of deformation. In Fig. 3, we show
the fineness and surface velocity fields for $m = 1.2$, $\hat{k}_A = 2.5$ and
$f_A = (0.05-0.75)10^{-4}$. At $f_A = 0.05 \times 10^{-4}$ the fineness field
shows a small perturbation above the baseline value (of 1) and the
associated boundary width is relatively small (see Figs 4c and e).
As we increase $f_A$ the maximum values of fineness increase as
expected ($\delta_{d,c} \sim (f_A/\hat{k}_A)^{1/(p-m)}$) and there is a slight broadening
of $\delta_{d,c}$ (see Fig. 4c). The convergent boundary width $\delta_{d,c}$ calculated
at $f_A = 0.05 \times 10^{-4}$ is smaller than at larger $f_A$ in Fig. 3, and
we can easily see that as $f_A$ is reduced the variations in the fine-
ness field get smoothed out until they disappear at $f_A = 0$. The
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Variations in surface velocity and fineness for variations in $f_A$ with $\hat{k}_A = 2.5$ and $m = 1.2$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** The range of plate-like behaviours observed (+, with connecting dashed lines) and predicted (solid lines) from our study for the case of variations in fineness inducing damage characterized by $f_A$. (a) Plateness versus $f_A$, (b) Plate velocity versus $f_A$, (c) Convergent boundary width versus $f_A$, (d) Convergent boundary strain-rate versus $f_A$, (e) Divergent boundary width versus $f_A$ and (f) Divergent zone strain-rate versus $f_A$.</figcaption></figure>

divergent boundary widths in Fig. 4(e) display less systematic behaviour as a function of $f_A$ than the convergent boundary widths, and in Fig. 3 appear to not change significantly. The increase in strain-rate as $f_A$ is increased, therefore, comes about due to the increase in $u_L$ at approximately constant $\delta_{c,d}$ (Figs 4d and f). Since increasing the fraction of fineness inducing damage does nothing to increase the rigidity of the plate (and $\delta_{c,d} \approx$ constant) the plateness remains approximately unchanged. Increasing the fraction of deformational work that goes into the generation of fineness damage at constant $\hat{k}_A$ and $m$, therefore, does not seem to improve the plate characteristics in any dramatic fashion contrary to previous consid-

<figure><img src="image_6.png" /><figcaption>**Figure 5.** The above plot shows the behaviours in boundary widths for variations in $f_A$ with the healing rate $\hat{k}_A$ kept constant for viscosity exponent $m = 2$ and two different healing rates. The grey inset is for $\hat{k}_A = 250$ and the other plot is for $\hat{k}_A = 2.5 \times 10^5$, but the axes are the same for both plots. For both healing rates we find that at small values of $f_A$ (where small depends upon healing rate) the system is in a stagnant lid regime. Further increases in $f_A$ forces the system to an episodic regime followed by a region characterized by stable plate behaviour. Increases in $f_A$ beyond the steady state regions results in numerically unresolvable phenomena.</figcaption></figure>

erations (Auth *et al.* 2003; Bercovici & Ricard 2005). However, in Bercovici & Ricard (2005) they looked at variations in $f_A/\sigma$ (with the surface tension $\sigma$ assumed small) and were possibly able to allow for more significant localization, though they primarily looked at generation of toroidal motion without characterizing boundary width variations. This study suggests that variations in the amount of fineness generating damage available in a convecting system at steady state produces the same amount of localization. Fig. 5 displays the results for variations in boundary widths for variations in $f_A$ at two different healing rates. For a given healing rate at small values of $f_A$ the system is in a stagnant lid regime, and increases in $f_A$ cause a transition to time dependent behaviour characterized by episodic plate-like behaviour. Upon further increasing $f_A$ the system enters a stable plate-like regime. For a given healing rate, there is little variation in localization in the stable plate regime. With that being said, increasing the fineness inducing damage characterized by $f_A$ does allow for a viscosity drop that increases plate mobility without decreasing boundary widths (and hence increasing localization).

### 3.2.2 Variations in healing rate $\hat{k}_A$

We now consider how variations in healing rate $\hat{k}_A$ affect the behaviours of convergent and divergent boundary widths (Figs 6c and e). Our results show that increasing the healing rate (at constant damage fraction $f_A$ and viscosity exponent $m$) significantly reduces plate velocity (Fig. 6b), has little effect on plateness (Fig. 6a), and has varied effects on convergent and divergent boundary width (Figs 6c and e). Increasing the healing rate leads to an increased plate viscosity which acts to decrease the plate velocity and boundary strain-rate. Convergent boundary widths show unsystematic behaviour to variations in healing rate except when $\hat{k}_A$ becomes small and $\delta_c$ generally decreases mildly. The trends of $\delta_d$ for all $m$ values suggest that the effect of increasing $\hat{k}_A$ initially leads to a reduced boundary width, but when $\hat{k}_A$ reaches a critical value divergent boundary width begins to increase and eventually reach an asymptote. The initial decrease in $\delta_d$ is consistent with our previous analysis that a decrease in the healing length scale $u_{max}/\hat{k}_A$ (dominated by the decrease in $u_{max}$) will lead to a reduction in $\delta_d$. Therefore, we see that an increase in healing can actually play a role in localization of divergent boundaries by causing the reduced
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The range of plate-like behaviours observed (+, with connecting dashed lines) and predicted (solid lines) from our study for the case of variations in healing rate characterized by $\hat{k}_A$ with $f_A = 0.5 \times 10^{-4}$. (a) Plateness versus $\hat{k}_A$, (b) Plate velocity versus $\hat{k}_A$, (c) Convergent boundary width versus $\hat{k}_A$, (d) Convergent boundary strain-rate versus $\hat{k}_A$, (e) Divergent boundary width versus $\hat{k}_A$ and (f) Divergent zone strain-rate versus $\hat{k}_A$.</figcaption></figure>

grain size at a ridge to substantially heal within a short distance from the ridge. When the healing rate reaches the critical value, the velocity (and boundary strain-rate) approach their small-magnitude asymptotes. At this point the lithosphere begins to approach the properties of a stagnant lid, and hence the boundary widths increase (divergent) or remain the same (convergent). Our expectation that increasing healing rate would lead to a moderation of localization while allowing for an increase in plate strength seems to be partially born out. While increasing healing rate does strengthen the plate interiors and allows for steady-state convection, it will eventually choke off surface motion and cause a rigid lid. The transition in convective behaviours would be similar in na-

ture to Fig. 5, except that stagnant lid (stable plate-like) behaviour would be at large (small) $k_A$ instead of small (large) $f_A$. An intriguing and surprising result is that increasing healing rate (in a specific range) may allow for localization at divergent plate boundaries via a reduction in the healing length scale.

### 3.2.3 Variations in $f_A$ and $\hat{k}_A$ with their ratio constant

We further consider the physics of advection in plate boundary development by maintaining fixed relative sizes of damage and healing. The previous sections showed that increasing $f_A$ resulted in an increase in plate velocity and boundary strain-rate, while increasing the healing rate ($\hat{k}_A$) resulted in a decrease in plate velocity and boundary strain-rate. Our results for constant $f_A/\hat{k}_A$ at various values of $m$ (note that the ratio $f_A/\hat{k}_A$ is different for each $m$ value, and is chosen to be the largest ratio that is numerically resolvable (e.g. $f_A/\hat{k}_A$ for $m = 1.2$ is $0.25 \times 10^{-4}$ while $f_A/\hat{k}_A$ for $m = 2$ is $0.002 \times 10^{-4}$)) are shown in Figs 7 and 8. We see that plate velocity remains fairly unchanged with variations in $f_A$ and $\hat{k}_A$ while $f_A/\hat{k}_A$ is kept constant (Fig. 8b), whereas we say $u_{max}$ grows with an increasing ratio of $f_A$ to $\hat{k}_A$ in Figs 4 and 6. Our results show that divergent boundary width undergoes a decline as $f_A$ and $\hat{k}_A$ increase before eventually asymptoting to its minimum value (Fig. 8e). This result again shows the phenomena discussed in the previous section where a decrease in the healing length scale ($u_{max}/\hat{k}_A$) may allow for localization of divergent boundaries. Contrary to the previous section though the plate velocity remains constant as $\hat{k}_A$ goes up and the model does not approach the stagnant lid case; therefore, increasing $f_A$ simultaneously with the healing rate allows for divergent boundaries to undergo further localization than when $f_A$ remains small. The asymptotic value for $\delta_d$ is eventually controlled by the competition between healing and deformational work when both fineness inducing damage and healing rate dwarf the advective component. We have already seen that convergent boundaries respond somewhat negligibly to variations in $f_A$, but here we see that $\delta_d$ actually grows for increasing healing rate and fineness inducing damage fraction (Fig. 8c). It seems that the effect of increasing $\hat{k}_A$ acts to distribute deformation by increasing the plate strength at old lithosphere, while the value of $f_A$ does not contribute significantly to weakening of the plate at subduction zones. These results again highlight the conclusion that the combination of plate velocity and healing may allow for narrow divergent boundaries, but increases in fineness inducing damage characterized by $f_A$ still remain insufficient for generating narrow convergent boundaries. In Fig. 9 we

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Variations in velocity and fineness for different $m$, $f_A$ and $\hat{k}_A$ values, but keeping $f_A/\hat{k}_A$ = constant.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

$$m = 1.2, f_A \mu_R / \hat{k}_A = .25 \qquad m = 1.8, f_A \mu_R / \hat{k}_A = .005$$

$$m = 1.4, f_A \mu_R / \hat{k}_A = .05 \qquad m = 2, f_A \mu_R / \hat{k}_A = .002$$

**Figure 8.** The range of plate-like behaviours observed (+, with connecting dashed lines) and predicted (solid lines) from our study for the case of simultaneous variations in $f_A$ and $\hat{k}_A$. (a) Plateness versus $f_A$, (b) Plate velocity versus $f_A$, (c) Convergent boundary width versus $f_A$, (d) Convergent boundary strain-rate versus $f_A$, (e) Divergent boundary width versus $f_A$ and (f) Divergent zone strain-rate versus $f_A$.

<figure><img src="image_10.png" /></figure>

**Figure 9.** The variation in plate boundary width as a function of damage fraction for $5 \times 10^{-5} \le f_A \le 0.1$. This case is for $m = 2$ and constant ratio of $f_A$ to $\hat{k}_A$ equal to $0.002 \times 10^{-4}$. We can see that the variations in plate boundary width saturate well before $f_A$ approaches 0.1.

see the effect of increasing $f_A$ to significantly larger values (and correspondingly $\hat{k}_A$) demonstrates that boundary widths still remain saturated. Therefore, our analysis of behaviours at smaller $f_A$ displays the variations brought about maintaining constant advection by fixing the ratio of damage to healing.

### 3.2.4 Variations in m

We finally consider how variations in viscosity exponent $m$ affect the plate-like behaviour of our model at constant values of $f_A$ and $\hat{k}_A$. The results shown in Fig. 10 are for various combinations of $f_A$ and $\hat{k}_A$, but all show that an increase in $m$ results in a decrease in boundary width for both $\delta_d$ and $\delta_c$ (Figs 10c and e), which is expected from our scaling analysis and previous observations (Auth et al. 2003; Bercovici & Ricard 2005). Our calculations also show that the plateness of the flow increases approximately linearly with

<figure><img src="image_11.png" /></figure>

$$f_A = 0.5 \times 10^{-4}, \quad \hat{k}_A = 25, 50, 100, 200, 250$$

**Figure 10.** The range of plate-like behaviours observed (+, with connecting dashed lines) and predicted (solid lines) from our study for the case of variations in viscosity exponent $m$ with $f_A = 0.5 \times 10^{-4}$. (a) Plateness versus $m$, (b) Plate velocity versus $m$, (c) Convergent boundary width versus $m$, (d) Convergent boundary strain-rate versus $m$, (e) Divergent boundary width versus $m$, and (f) Divergent zone strain-rate versus $m$.

$m$ (Fig. 10a), irrespective of the ratio $f_A$ to $\hat{k}_A$. Increases in the viscosity exponent, therefore, also increase the rigidity of the plate while also decreasing the boundary widths. As $m$ is increased, the boundary zone strain-rate increases dramatically (Figs 10d and f), primarily due to the decrease in boundary width, but strain-rate is also influenced by the corresponding increase in plate velocity. Allowing for greater viscosity sensitivity to grain size by increasing $m$, therefore, allows for significantly greater plate-like behaviours by increasing plate velocity and plateness, while narrowing both divergent and convergent plate boundaries.

In Fig. 11, we show the relevant plate characteristics for a case where $m = 2.5$, $f_A = 0.5 \times 10^{-4}$, and $\hat{k}_A = 750$ at steady state along with its time series. The time series values of plateness, Nusselt number and maximum fineness show the transition from our initial state, through the period of transient oscillations to an eventual steady state. The velocity field is very plate-like, with narrowly deforming boundaries and a rigid plate interior. The value of fineness at the subduction zone is approximately five, while the minimum value of fineness in the plate interior is about 0.05 for a range of fineness of about two orders of magnitude. Given our reference grain size of 1 mm, this would correspond to a range of grain sizes of about 200 µm to 20 mm. While this range of values for grain size may be possible within the Earth, we would expect the grain size in a shear zone and plate interior to be somewhat reduced compared to our calculated values (Evans et al. 2001). In our efforts to find damage parameters that are able to produce plate like flows though it is encouraging that successful results don't require unreasonable values of material parameters, such as grain size.

## 3.3 Combined void- and fineness-generating damage

In this final section we consider the combined effects of both void- and fineness-generating damage mechanisms on plate-like
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Results for parameter set $m = 2.5$, $f_A = 0.5 \times 10^{-4}$ and $\hat{k}_A = 750$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Results for case when both void and fineness generating damage are generated by deformational work. In these plots $\hat{k}_A = 2.5 \times 10^6$, $m = 2$ and $\gamma = 10^{-2}$.</figcaption></figure>

behaviour in the model. We have so far found that void generating is insufficient by itself to produce plate like flow in the current two-phase framework, which is in line with previous work (Bercovici & Ricard 2005). On the other hand, fineness generating damage has shown to be successful at producing plate like flow under certain circumstances in this study as well as previously (Bercovici & Ricard 2005). Fig. 12 shows results for a case with both void and fineness generating damage, the void- and fineness-generating end member cases are presented alongside for comparison. Previously we found that employing a damage rheology with only void generating damage and $\gamma = 10^{-2}$ produced very unplate-like flows (see Fig. 12). However, when these same parameters are employed with moderate fineness generating damage ($f_\phi = f_A = 0.25$, $\hat{k}_A = 2.5 \times 10^6$, $m = 2$) the result is localized plate boundaries and significantly more plate-like flow. The combined case in fact becomes too localized to numerically resolve

subsequent to frames shown. The reason for the ultra localization in this case is due to the $f_A/\eta(\phi)$ term in eq. (15). The void-generating damage allows for the porosity field to undergo larger variations than when $f_\phi = 0$, and hence $\eta(\phi)$ is reduced in the convergent and divergent zones. The reduction in $\eta(\phi)$ results in an increase in the effective amount of fineness generating damage, $f_A/\eta(\phi)$. Physically, the reduction in $\eta(\phi)$ (or interfacial area density) allows for the continued input of damage to have a more significant impact than if the system was already heavily damaged, since inputting deformational work into a material with little damage results in a greater viscosity drop than if the material was already significantly deformed. While our results for boundary width variations in $f_A$ (with $f_\phi = 0$ at constant $\hat{k}_A$ and $m$) (Fig. 4) didn't show any significant localization it seems that a larger effective $f_A$ (with $f_\phi > 0$) is capable of generating much greater plate-like flow. When $\gamma$ is increased to $10^3$ (with the same fineness parameters)
<figure><img src="image_14.png" /><figcaption>**Figure 13.** The effect of combined void and fineness generating damage on plate-like behaviour for different values of total damage and fineness healing. For this plot $m = 2$, $f = 0.05$, $0.5$, $\hat{k}_A = 2.5 \times 10^5$, $2.5 \times 10^6$ and $\gamma = 10^{-2}$. The quantity $f$ is defined to be the sum of $f_A$ and $f_\phi$.</figcaption></figure>

we eventually force the system to be similar to the case when $f_\phi = 0$.

The results of allowing combined void- and fineness-generating damage in our convective model is summarized in Fig. 13. For some cases with combined forms of damage the system became numerically unresolvable (i.e. solution became too localized before reaching steady state), hence we choose to measure the maximum strain-rate when the convergent boundary width was less than three gridpoints. For each choice of $m$, $k_A$ and $\gamma$ the plate velocity is generally the same, so increasing strain-rate $[(u/\delta)_{\text{max}}]$ is indicative of a reduction in boundary width and hence localization. The results show that an increase in plate-like behaviour, specifically localization, is achieved by allowing for both forms of damage when the fineness generating parameters by themselves are associated with stable plate-like behaviour ($\hat{k}_A = 2.5 \times 10^5$, $f = 0.05$ and $\hat{k}_A = 2.5 \times 10^6$, $f = 0.5$ in Fig. 13). In contrast, the combination of void and fineness generating damage, when the fineness parameters are associated with episodic or stagnant lid behaviour ($\hat{k}_A = 2.5 \times 10^6$, $f = 0.05$ in Fig. 13), does not result in an increase in plate localization or mobilization. Combining the two types of damage mechanisms can result in more significant localization, but this seems to occur only when the fineness parameters permit stable plate-like solutions. Addition of void generating damage does not enhance plate like behaviour unless the fineness parameters allow for it. As previously discussed, increasing the parameter $\gamma$ acts to force the system back to the case when $f_\phi = 0$ and decreasing $\gamma$ causes the convective system to localize at a faster rate. These results suggest that an interesting interplay between void and fineness generating damage may exist that could facilitate shear localization within the Earth. From a mineral physics perspective microcracks distributed in a poly-crystalline material may act to pin grain boundaries, therefore, inhibiting grain growth. If grain growth is stunted due to the presence of pores (by $f_\phi \neq 0$) then grain size will remain smaller and viscosity will correspondingly remain small. Results from looking at natural samples has shown that grain size reduction can lead to the formation of brittle faults and shear localization (Jin et al. 1998). Brittle failure may lead to significant grain size reduction as well, though the reduction process may be more related to generation of fault gouge than grain size reduction via dynamic recrystallization or similar processes. The interaction of brittle and ductile deformation mechanisms remains a poorly understood rheological phenomena from both an experimental and theoretical perspective (Kohlstedt et al. 1995).

# 4 DISCUSSION AND CONCLUSION

## 4.1 Earth-like parameters and our model

Plate generation necessitates breaking the strongest part of the lithosphere; we would therefore, like to know if the range of parameters in our thin-sheet model is associated with the strongest portion of the lithosphere [i.e. brittle–ductile transition (Kohlstedt et al. 1995)]. One way of estimating this is through the grain growth rate $k_A$ in eq. (15). From mineral physics we know that the grain growth rate in olivine and other minerals is highly temperature dependent (Karato 1989). The growth rate has an Arrhenius type functional form that can be solved for temperature such that

$$k_A = k_0 \exp \left( \frac{-H}{RT} \right) \Rightarrow T = \frac{H}{R \ln \left( \frac{k_0}{k_A} \right)}, \quad (45)$$

where $k_0$ is the reference grain growth rate, $H$ is the activation enthalpy, and $R$ is the gas constant. Given the values of $\hat{k}_A$ we determined to be capable of producing plate-like flows we can then get an estimate of the lithospheric temperature controlling the rheology in our model. Using values of $k_0 = 10^{-8} \text{ m}^2 \text{ s}^{-1}$ and $H = 2 \times 10^5 \text{ J mol} \text{ K}^{-1}$ (Karato 1989) and $\hat{k}_A = 250–1000(2.3–9.0 \times 10^{-21} \text{ m}^2 \text{ s}^{-1}$ in dimensional form) for $m = 2$ and $f_A \mu_R = 0.5$ from our calculations we arrive at a temperature range of $T = 800–900 \text{ K}$. This range of values for temperature will put us in the brittle–ductile regime for most of the oceanic lithosphere. Our thin-sheet damage rheology is, therefore, associated with the strongest region of the lithosphere, which is what is desired in a thin-sheet model. The kinetics of grain growth at such low temperatures might be significantly different than eq. (45) predicts, but experimental data on single-phase grain growth likely predict much faster grain growth than actually occurs in the mantle due to the inhibition of grain growth by secondary phases (Evans et al. 2001). Therefore, our temperature estimate might be a lower limit of the temperature or depth represented by our vertically averaged rheology.

## 4.2 Final thoughts

The goal of this study is to test our hypothesis that two-phase damage theory will allow for shear localization leading to the development of plate tectonics on Earth. We have employed a 2-D model of mantle convection in order to test the plate generating capabilities of two-phase damage theory. The primary characteristics of plate tectonics on Earth that we have attempted to address in this study include the formation of narrow boundaries at both convergent and divergent zones, rigid plate interiors, and temporally persistent plate-like velocity profiles. We have attempted to understand how the different factors (i.e. advection, deformational work, and healing) affect the rheological variables (fineness and porosity) to produce the aforementioned plate-like features.

We have explored how each form of damage, void and fineness generating (or grain size reducing), by itself and in tandem, can facilitate plate-like behaviour in a convectively driven system. Our results for void and fineness generating damage by themselves have shown some similar consequences to separate problems as previous work (Bercovici & Ricard 2005), namely that void generating damage was unsuccessful while fineness generating damage was quite successful at inducing plateness in this study and toroidal motion in Bercovici & Ricard (2005). An important part of this study has been to understand how the range of different parameters that govern grain size (fineness) evolution affect the plate
characteristics we find. We have found that increasing fineness in-
ducing damage $f_A$ (for fixed healing rate $k_A$ and viscosity expo-
nent $m$) does not do much to improve the plate characteristics in
our model, though allowing for void-generating damage in tandem
with fineness-generating damage dramatically changes this result.
Increases in the healing term $k_A$ can lead to an important reduction
in divergent boundary width due to the reduction in the healing
length scale. This behaviour may also be important for the forma-
tion of narrow rifts in that plate thinning under tension may allow
for significant healing outside of the failure region thereby forming
a localized rift zone. Increases in the viscosity exponent $m$ lead to
the formation of narrow plate boundaries (both convergent and di-
vergent) as well as increases in the rigidity of the plate interior, and
the values of $m$ explored fall within the range of values expected for
mantle rheologies (Karato & Wu 1993). These results suggest that
the mechanism of localization at convergent and divergent bound-
aries need not originate from the same process. Our model does not
incorporate the generation of melts which will strongly influence
the strength of ridges as well (Tackley 2000a,b), and melt forma-
tion obviously leads to different localization mechanism than what
affects subduction zones. While grain size sensitive creep (coupled
with crack formation) may be needed to form narrow convergent
boundaries, the generation of narrow divergent boundaries may
only need some small damage input coupled with strong healing
near the ridge. Allowing for both forms of damage causes a feed-
back mechanism which allows for the most significant localization
of plate boundaries in our study. The feedback mechanism begins
with void-generation allowing for greater variability in porosity, the
range of $\eta(\phi)$ (or interfacial area density) increases, which leads
to greater damage input into increasing fineness (hence reduction
in boundary viscosity), which focuses the variations in the poros-
ity field, hence completing the cycle that facilitates the formation
of narrow plate boundaries. While many studies have considered
the effect of cracking and grain size reduction on plate generation
or convection separately (Hall & Parmentier 2003; Bercovici &
Ricard 2005), our study suggests that it is the interaction of these
two forms of damage which allows for the formation of plate tec-
tonics on Earth. Further investigation into the coupling between
void- and fineness-generating damage in a depth-dependent prob-
lem, hence allowing us to follow the plate boundary formation
process in a vertical section of the lithosphere, is a necessary next
step.

While our model does include the zeroth order effect of temper-
ature on the viscosity of the lithosphere (i.e. lithospheric viscosity
is much greater than Newtonian mantle viscosity), the effect of a
laterally varying temperature on the lithospheric dynamics is an im-
portant consideration to look into for the future. A laterally variable
temperature in the plate would likely allow for the formation of even
narrower divergent plate boundaries. The area of ridge formation is
significantly warmer than assumed in the model, and including this
effect would lead to a lower viscosity as well as increase the amount
of healing in the grain size evolution equation. The viscosity of
the region of old lithosphere would be similar to the zeroth order
viscosity assumed in our model and, therefore, would likely be little
affected by allowing for lateral temperature variations.

While a thin-sheet formulation for lithospheric dynamics cap-
tures a significant component of lithospheric behaviour, allowing
for depth-dependent rheologies might allow for a wider range of
plate generating scenarios. As we noted previously the likelihood
of void-generating damage manifesting itself at mid to deep litho-
spheric depths is small, but it is entirely possible that localization in
the lithosphere may vary as a function of depth. Allowing for depth-

dependent rheologies would enable us to explore the possibilities of
fault formation and progression in a more complete manner. For ex-
ample, localization may be initiated by grain size reduction at mid to
deep lithospheric depths leading to failure at shallower depths via
void generation completing the formation of a lithospheric fault.
These considerations are an important step for future analyses, but
the current thin-sheet model does allow for a fairly comprehensive
understanding of how two-phase damage theory could facilitate the
formation of tectonic plates on Earth.

## ACKNOWLEDGMENTS

The authors benefited from discussions with Shun-ichiro Karato and
Garrett Leahy, as well as the feedback in the Yale mantle/core semi-
nar. The suggestions from two anonymous reviewers also helped to
clarify many different aspects of this paper. Support was provided
by the National Science Foundation (NSF, grant EAR-0105269) and
the Centre National de la Recherche Scientifique (CNRS).

## REFERENCES

Austin, N. & Evans, B., 2007. Paleowattmeters: a scaling relation for dy-
namically recrystallized grain size, *Geology*, **35**(4), 343–346.

Auth, C., Bercovici, D. & Christensen, U., 2003. Two-dimensional convection with a self-lubricating, simple-damage rheology, *Geophys. J. Int.*, **154**, 783–800.

Bercovici, D., 1993. A simple model of plate generation from mantle flow, *Geophys. J. Int.*, **114**, 635–650.

Bercovici, D., 1995. A source-sink model of the generation of plate tectonics from non-Newtonian mantle flow, *J. geophys. Res.*, **100**, 2013–2030.

Bercovici, D., 1998. Generation of plate tectonics from lithosphere–mantle flow and void-volatile self-lubrication, *Earth planet. Sci. Lett.*, **154**, 139–151.

Bercovici, D., 2003. The generation of plate tectonics from mantle convection, *Earth planet. Sci. Lett.*, **205**, 107–121.

Bercovici, D. & Ricard, Y., 2003. Energetics of a two-phase model of lithospheric damage, shear localization and plate-boundary formation, *Geo-phys. J. Int.*, **152**, 581–596.

Bercovici, D. & Ricard, Y., 2005. Tectonic plate generation and two-phase damage: void growth versus grain size reduction, *J. geophys. Res.*, **110**, 1–18.

Bercovici, D., Ricard, Y. & Richards, M., 2000. The relation between mantle dynamics and plate tectonics: a primer, in *History and Dynamics of Global Plate Motions, Geophys. Monogr. Ser.*, Vol. 121, pp. 5–46, eds Richards, M.A., Gordon, R. & van der Hilst, R., Am. Geophys. Union, Washington, DC.

Bercovici, D., Ricard, Y. & Schubert, G., 2001a. A two-phase model of compaction and damage, 1. General theory, *J. geophys. Res.*, **106**(B5), 8887–8906.

Bercovici, D., Ricard, Y. & Schubert, G., 2001b. A two-phase model of compaction and damage, 3. Applications to shear localization and plate boundary formation, *J. geophys. Res.*, **106**(B5), 8925–8940.

Cawood, P.A., Kroner, A. & Pisarevsky, S., 2006. Precambrian plate tectonics: criteria and evidence, *GSA Today*, **16**(7), 4–11.

Christensen, U. & Harder, H., 1991. Three-dimensional convection with variable viscosity, *Geophys. J. Int.*, **104**, 213–226.

Dumoulin, C., Bercovici, D. & Wessel, P., 1998. A continuous plate-tectonic model using geophysical data to estimate plate margin widths, with a seismicity based example, *Geophys. J. Int.*, **133**, 379–389.

England, P. & McKenzie, D., 1982. A thin viscous sheet model for continental deformation, *Geophys. J. R. astr. Soc.*, **70**, 295–321.

Evans, B., Renner, J. & Hirth, G., 2001. A few remarks on the kinetics of static grain growth in rocks, *Int. J. Earth Sci.*, **90**, 88–103.
Gurnis, M., Zhong, S. & Toth, J., 2000. On the competing roles of fault reactivation and brittle failure in generating plate tectonics from mantle convection, in *History and Dynamics of Global Plate Motions, Geophys. Monogr. Ser.*, Vol. 121, pp. 73–94, eds Richards, M.A., Gordon, R. & van der Hilst, R., Am. Geophys. Union, Washington, DC.

Hager, B. & O'Connell, R., 1981. A simple global model of plate dynamics and mantle convection, *J. geophys. Res.*, **86**, 4843–4867.

Hall, C. & Parmentier, E., 2003. Influence of grain size evolution on convective stability, *Geochem. Geophys. Geosystems* (G3), 1, 2003GC000308.

Jin, D., Karato, S. & Obata, M., 1998. Mechanisms of shear localization in the continental lithosphere: Inference from the deformation microstructures of peridotites from the Ivrea zone, northwestern Italy, *J. Struct. Geol.*, **20**, 195–209.

Karato, S., 1983. Grain-size distribution and rheology of the upper mantle, *Tectonophysics*, **104**, 155–176.

Karato, S., 1989. Grain growth kinetics in olivine aggregates, *Tectonophysics*, **168**, 255–273.

Karato, S. & Wu, P., 1993. Rheology of the upper mantle: a synthesis, *Science*, **260**, 771–778.

Karato, S., Toriumi, M. & Fujii, T., 1980. Dynamic recrystallization of olivine single crystals during high temperature creep, *Geophys. Res. Lett.*, **7**, 649–652.

Katz, R., Spiegelman, M. & Holtzman, B., 2006. The dynamics of melt and shear localization in partially molten aggregates, *Nature*, **442**, 676–679.

King, S., 2001. Subduction zones: observations and geodynamic models, *Phys. Earth planet. Int.*, **127**, 9–24.

Kohlstedt, D., Evans, B. & Mackwell, S., 1995. Strength of the lithosphere: Constraints imposed by laboratory experiments, *J. geophys. Res.*, **100**, 17587–17602.

Lemery, C., Ricard, Y. & Sommeria, J., 2000. A model for the emergence of thermal plumes in Rayleigh-Benard convection at infinite Prandtl number, *J. Fluid Mech.*, **414**, 225–250.

McKenzie, D., 1984. The generation and compaction of partially molten rock, *J. Petrol.*, **25**, 713–765.

Moresi, L. & Solomatov, V., 1998. Mantle convection with a brittle lithosphere: thoughts on the global tectonic style of the Earth and Venus, *Geophys. J. Int.*, **133**, 669–682.

Ogawa, M., 2003. Plate-like regime of a numerically modeled thermal convection in a fluid with temperature-, pressure-, and stress-history-dependent viscosity, *J. geophys. Res.*, **108**, 2067, doi:10.1029/2000JB000069.

Ribe, N., 1992. The dynamics of thin shells with variable viscosity and the origin of toroidal flow in the mantle, *Geophys. J. Int.*, **110**, 537–552.

Ricard, Y. & Bercovici, D., 2003. Two-phase damage theory and crustal rock failure: the theoretical ‘void’ limit, and the prediction of experimental data, *Geophys. J. Int.*, **155**, 1057–1064.

Ricard, Y. & Froidevaux, C., 1986. Stretching instabilities and lithospheric boudinage, *J. geophys. Res.*, **91**, 8314–8324.

Ricard, Y., Bercovici, D. & Schubert, G., 2001. A two-phase model of compaction and damage, 2, Applications to compaction, deformation, and the role of interfacial surface tension, *J. geophys. Res.*, **106**(B5), 8907–8924.

Schubert, G., Turcotte, D. & Olson, P., 2001. *Mantle Convection in the Earth and Planets*, Cambridge Univ. Press, Cambridge, UK.

Spiegelman, M., 1993a. Flow in deformable porous media, part 1, Simple analysis, *J. Fluid Mech.*, **247**, 17–38.

Spiegelman, M., 1993b. Physics of melt extraction: theory, implications and applications, *Phil. Trans. R. Soc. London, Ser. A*, **342**, 23–41.

Tackley, P., 1998. Self-consistent generation of tectonic plates in three-dimensional mantle convection, *Earth planet. Sci. Lett.*, **157**, 9–22.

Tackley, P., 2000a. Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations, 1. Pseudoplastic yielding, *Geochem. Geophys. Geosystems* (G3), **1**, 2000GC000036.

Tackley, P., 2000b. Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations, 2. Strain

weakening and asthenosphere, *Geochem. Geophys. Geosystems* (G3), 1, 2000GC000043.

Wdowinski, S., O'Connell, R. & England, P., 1989. A continuum model of continental deformation above subduction zones: Application to the Andes and the Aegean, *J. geophys. Res.*, **94**, 10331–10346.

Weinstein, S., 1996. Thermal convection in a cylindrical annulus with a non-Newtonian outer surface, *PAGEOPH*, **146**, 551–572.

Weinstein, S. & Olson, P., 1992. Thermal convection with non-Newtonian plates, *Geophys. J. Int.*, **111**, 515–530.

## APPENDIX: SCALING ANALYSIS: PREDICTIONS VERSUS RESULTS

Given the extensive numerical data generated, we test how well the scaling analyses developed in the previous sections fit the data. In Fig. A1, we show the results of plotting $A_o$ (the maximum steady state value of fineness, see eq. 37) at both ridges and subduction zones versus eq. (43) for each numerical experiment. The fit of the data to the prediction is quite good, and indicates that our assumption that stress (deformational work) is constant in a given environment (i.e. divergent vs. convergent) is reasonably accurate. The slopes of the predicted fineness magnitudes are actually different for ridges and subduction zones owing to the fact that the second invariant of stress between these two different regions is different. These results suggest that the steady-state amplitude of fineness at subduction zones and ridges is governed by the competition between deformational work and healing, and advection plays no role in determining the amplitude.

The results of comparing our numerical data for $\delta_c/u_{max}$ and $\delta_d/u_{max}$ with our scaling analysis is presented in Fig. A2(a). In Fig. A2(b) we plot ratio of the numerical result for $\delta_{c,d}/u_{max}$ to the predicted value for each point in parameter space, since we expect

<figure><img src="image_15.png" /><figcaption>**Figure 14.** The values of fineness plotted on the ordinate come from the numerical experiments, while the values of fineness plotted on the abscissa are from the scaling analyses eq. (43). The results from the numerical experiments and the scaling analyses are for the maximum values of steady-state fineness ($A_o$) in both subduction zones and ridges.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** (a) Ordinate values are boundary widths over $u_{max}$ (or one over the strain-rate) that we determined from numerical experiments plotted against the values on the abscissa which are from the scaling analyses. (b) Ordinate values are the ratio of inverse strain-rate ($\delta_{c,d}/u_{max}$) for the scaling analysis to the numerical experiments plotted against the run number. This graph displays the constant scale that differs between the scaling predictions and numerical experiments.</figcaption></figure>
our scaling analysis to only fit the numerical data up to some multi-
plicative constant (e.g. a geometrical factor of order 1 not included
in a scaling analysis). The fit for convergent boundaries is especially
good, but both boundaries exhibit the general features suggested by
the scaling in eq. (40) as demonstrated by the linear relationship

between the numerical experiments and analytic prediction (where
the analytic prediction has been multiplied by the constant deter-
mined in Fig. A2b). The regression line fits all of the data points to
within approximately 15 per cent which we will assume allows us
to assign some validity to the analytic expression we derived.