2. Model

The model envisions an axisymmetric viscous shell with density $\rho _1=\rho +\delta \rho$ and viscosity $\eta _1$ immersed in a spherical viscous planet with density $\rho _0=\rho$ and viscosity $\eta _0$ (figure 1). The flow in both fluids is governed by the incompressible Stokes equations,

(2.1a,b)\begin{equation} -\boldsymbol{\nabla}p+\eta_i\nabla^2\boldsymbol{u}+\rho_i\boldsymbol{g}= \boldsymbol{0},\quad\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{equation}

where $p$ is the pressure, $\boldsymbol {u}$ the velocity and $\rho _i\boldsymbol {g}$ is the gravitational force per unit volume. The gravitational acceleration is directed radially toward the centre of the planet, and its magnitude $g$ is treated as a constant for simplicity ($g$ is in any case nearly constant throughout the Earth's mantle). The outer surface $r = R_0$ of the planet is free slip, i.e. the normal component of velocity and the shear traction both vanish there. Our model lacks the effectively inviscid core of radius $\approx 0.54 R_0$ that exists in the Earth. This neglect is deliberate: it removes a non-essential length scale from the problem that would complicate physical interpretation if it were retained. Neglecting the core also makes it possible to derive analytical Green functions for the flow due to a point force, which in turn permit us to employ an efficient form of the BEM to solve our problem. We also neglect the planet's deviation from perfect sphericity due to the centrifugal force associated with its rotation, which is typically small compared with $R_0$.

Figure 1.

Definition sketch of the model. Right portion (not to scale): an axisymmetric viscous shell with thickness $h$, density $\rho + \delta \rho$ and viscosity $\eta _1$ is immersed in a spherical viscous planet with radius $R_0$, density $\rho$ and viscosity $\eta _0$. The outer surface of the planet is impermeable and free of shear traction. The shell comprises a spherical ‘plate’ $\theta \leqslant \theta _t$ and a downward-dipping ‘slab’ $\theta _t\leqslant \theta \leqslant \theta _s$. The lubrication layer above the plate has thickness $d$, and the radius of the unsubducted part of the plate's midsurface (dotted line) is $R\equiv R_0 - h/2 - d$. The gravitational vector $\boldsymbol {g}$ is directed radially toward the centre of the planet. Left portion: the same shell drawn to a scale appropriate for the Earth, with $h = 100$ km, $d=20$ km and $R_0 = 6370$ km.

The first step in the BEM approach is to specify the initial shape of the shell. As indicated in figure 1, the model shell comprises two pieces that join together smoothly at a colatitude $\theta = \theta _t$. The inner piece $\theta \leqslant \theta _t$ is a circular spherical cap of constant thickness $h$, separated from the free-slip surface by a thin ‘lubrication layer’ of thickness $d$. The radius of the midsurface of the cap is therefore $R = R_0 - d - h/2$. The outer piece $\theta _t\leqslant \theta \leqslant \theta _s$ is a downward-dipping ring of arcwise length $l$ whose midsurface (dotted line in figure 1) is located at the radius

(2.2a,b)\begin{equation} r(\theta) = R (1 - b \zeta^3 - c\zeta^4), \quad \zeta = \frac{\theta - \theta_t}{\theta_s - \theta_t} , \end{equation}

where $b$ and $c$ are constants. The thickness of the ring is $h$ everywhere outside the rounded rim just beyond $\theta = \theta _s$. Following geophysical parlance, we shall call the cap $\theta < \theta _t$ the ‘plate’, the circle $\theta = \theta _t$ the ‘trench’ and the downward-dipping ring $\theta _t\leqslant \theta \leqslant \theta _s$ the ‘slab’ (these terms explain the choice of subscripts on $\theta _t$ and $\theta _s$).

The mathematical form of (2.2a,b) ensures that the local slope and local curvature (in the $\theta$-direction) of the slab match those of the plate at the trench $\zeta = 0$. The values of $b$ and $c$ are determined by imposing two additional constraints on (2.2a,b). The first is that the dip $\varphi$ of the leading end of the slab's midsurface relative to the local horizontal is a specified value $\varphi _s$. The second constraint is that the curvature at the leading end of the midsurface is $-1/R$. The reason for imposing this constraint is that the slab's leading edge is free, i.e. the bending moment is zero there. Because the bending moment in a viscous shell is proportional to the rate of change of curvature of the midsurface, the curvature at the leading edge of the slab should not deviate from the value $-1/R$ for an undeformed spherical shell. Explicitly, the foregoing two constraints are

(2.3)\begin{gather} \left.\frac{r}{\left(r^2 + (\partial_{\theta} r\right)^2)^{1/2}}\right |_{\zeta = 1} = \cos\varphi_s, \end{gather}

(2.4)\begin{gather}\left. -\frac{r^2 + 2 (\partial_{\theta} r)^2 - r\partial^2_{\theta\theta} r} {\left(r^2 + (\partial_{\theta} r\right)^2)^{3/2}}\right |_{\zeta = 1} ={-}\frac{1}{R}. \end{gather}

Equations (2.3) and (2.4) are a rather complicated quartic system of algebraic equations for $b$ and $c$, which we solved numerically for given values of $R$, $\varphi _s$ and $\theta _s - \theta _t$.

A few words of explanation are in order concerning the lubrication layer above the plate in figure 1, which corresponds to the ‘sticky air’ layer mentioned earlier. The function of this layer is twofold: it creates large normal stresses that support the plate and prevent it from sinking; and it applies a negligible shear stress to the upper surface of the plate. The latter surface is thus effectively a free surface whose radial velocity is not constrained, allowing the plate to bend in a realistic way. The advantage of the lubrication layer formulation is that it obviates the need to deal with the numerically troublesome complexities of a true free surface, including the triple point at the trench where the shell, the ambient fluid and the overlying ocean would meet. Li & Ribe (Reference Li and Ribe2012) have shown that the predictions of Cartesian three-dimensional subduction models with a lubrication layer agree quantitatively (not just qualitatively) with laboratory experiments in which the plate is prevented from sinking by surface tension. This shows that two different mechanisms for preventing the sinking of the plate lead to identical predictions, inspiring confidence in the lubrication layer approach.

A final length scale that we shall need is the ‘bending length’ $l_b$, which plays a crucial role in the subsequent scaling analysis. When the Stokes equations are solved for an immersed shell like that shown in figure 1, the resulting velocity field will be non-zero everywhere inside the shell. Because the shell is thin, however, that velocity field can be described as a combination of bending and stretching. Based on previous experience with two-dimensional shells (Ribe Reference Ribe2010), we anticipate that the rate of deformation will be dominated by bending throughout the slab and in an adjoining portion of the plate. We call the arcwise length of this part of the shell the ‘bending length’ $l_b$. In geophysical terms, the bending length is the sum of the slab length $l$ and the length $l_b - l$ of the near-trench portion of the plate where significant upward motion (‘flexural bulging’) occurs. The bending length is indicated schematically in the right-hand portion of figure 1. Because the bending length is a dynamic length scale that arises physically from solving the Stokes equations, it is of a fundamentally different character than the geometric length scales $h$, $d$ and $l$.

3. Dimensional and scaling analysis

3.1. Dimensionless groups

We begin by performing a dimensional analysis to get an idea of the ‘size’ of the problem we face. The first task is to choose a characteristic output parameter of the model to serve as the target of the analysis. Following our earlier studies of subduction in two- and three-dimensional Cartesian geometries (Ribe Reference Ribe2010; Li & Ribe Reference Li and Ribe2012), we choose the vertical (i.e. radial) component $V$ of the velocity of the leading edge of the slab. This choice is based on the fact that in Cartesian geometry, $V$ satisfies a simpler scaling law than other possible choices such as the norm or the horizontal component of the leading-edge velocity (Ribe Reference Ribe2010).

Referring to figure 1, we see that $V$ depends on the planet radius $R_0$, the shell thickness $h$, the lubrication layer thickness $d$, the slab length $l$, the slab dip $\varphi _s$, the angle $\theta _t$ subtended by the plate, the mantle viscosity $\eta _0$, the shell viscosity $\eta _1$ and the buoyancy $g\delta \rho$. The gravitational acceleration and the densities of the shell and the mantle appear only in the combination $g\delta \rho$ because all velocities in gravity-driven Stokes flow are linearly proportional to the driving buoyancy force.

The foregoing list comprises ten parameters, three of which have independent dimensions. Buckingham's $\varPi$-theorem (Buckingham Reference Buckingham1914) therefore tells us that seven independent dimensionless groups can be formed from our original ten parameters. Incorporating $V$ into only one of the groups for convenience, we see that the scaling law satisfied by $V$ must have the general form

(3.1)\begin{equation} \frac{V \eta_0}{hl g\delta\rho} = f_0 \left( \frac{h}{R_0}, \frac{h}{l}, \frac{d}{h}, \frac{\eta_1}{\eta_0}, \theta_t, \varphi_s \right) , \end{equation}

where $f_0$ is an unknown function. The seven dimensionless groups could of course have been defined differently, but the choices in (3.1) are as good as any.

While dimensional analysis has reduced the number of independent parameters substantially, we are still faced with an unknown function of six arguments. Characterising such a function completely using numerical simulations is not practical, and would not provide any insight even if it were. However, we have not yet exploited the physics of our problem, which we now do using a scaling analysis based on the theory of thin viscous shells.

3.2. Physical scalings

Our scaling analysis focuses on the three forces that act on the bending portion (of length $l_b$) of the shell. These are the buoyancy force $F_b$, the external viscous force $F_{ext}$ exerted on the shell by the surrounding fluid and the internal viscous force $F_{int}$ that resists bending of the shell. The goal is first to find scalings for each of these in terms of the problem parameters, and then use them to find the scaling of the sinking rate $V$ of the slab.

The buoyancy force acting on the slab (per unit circumferential length) scales as

(3.2)\begin{equation} F_b\sim l h g\delta\rho. \end{equation}

The buoyancy force is proportional to $l$ rather than $l_b$ because the buoyancy of the shell in the flexural bulging region is exactly compensated by normal stresses in the overlying lubrication layer (Ribe Reference Ribe2010).

The viscous traction (force per unit area) applied to the shell by the outer fluid is $\sigma _{ext} \sim \eta _0 V/l_b$, where $V > 0$ is the downward radial velocity at the tip of the slab. Integrating this over the bending length, we find the external force per unit circumferential length acting on the shell

(3.3)\begin{equation} F_{ext}\sim \eta_0V. \end{equation}

Now we turn to the internal viscous force $F_{int}$ that resists deformation of the shell. This force is just the integral $N_1$ of the radial shear stress acting on a cross-section of the shell located at the plateward end of the bending length. To estimate this force, we exploit the theory of thin elastic shells (Novozhilov Reference Novozhilov1959) together with the Stokes–Rayleigh analogy between incompressible elasticity and slow viscous flow (Rayleigh Reference Rayleigh1945). According to this analogy, elastic shell theory can be transformed into its viscous equivalent by interpreting displacements as velocities and making the transformations $E\to 3\eta$ and $\sigma \to 1/2$, where $E$ is Young's modulus, $\eta$ is the viscosity and $\sigma$ is Poisson's ratio. From (7.8) of Novozhilov (Reference Novozhilov1959), we have

(3.4)\begin{equation} F_{int} \equiv N_1 = \frac{1}{A_1 A_2} \left[ \partial_{\theta} (A_2 M_1) - M_2\partial_{\theta} A_2 \right], \end{equation}

where

(3.5a,b)\begin{equation} A_1 = \left[ r^2 + (\partial_{\theta} r)^2\right]^{1/2} \quad\mathrm{and}\quad A_2 = r\sin\theta \end{equation}

are the Lamé parameters of the deformed axisymmetric shell and $M_1$ and $M_2$ are bending moments. The explicit expressions for these are

(3.6a,b)\begin{equation} M_1 = \frac{\eta h^3}{3}\left( \kappa_1 + \frac{1}{2}\kappa_2 \right), \quad M_2 = \frac{\eta h^3}{3}\left( \kappa_2 + \frac{1}{2}\kappa_1 \right) , \end{equation}

where

(3.7a,b)\begin{equation} \kappa_1 ={-} \frac{1}{A_1}\partial_{\theta} \left( \frac{1}{A_1}\partial_{\theta} W + K_1 U \right), \quad \kappa_2 ={-} \frac{1}{A_1 A_2}\partial_{\theta} A_2 \left( \frac{1}{A_1}\partial_{\theta} W + K_1 U \right), \end{equation}

and $U$ and $W$ are the velocity components tangential to and normal to the midsurface, respectively. The quantity $K_1$ in (3.7a,b) is the curvature of a line of constant longitude; its explicit expression is the left-hand side of (2.4). Physically, the quantities $-\kappa _1$ and $-\kappa _2$ are the rates of change of curvature of the midsurface in the $\theta$- and $\phi$-directions, respectively, due to deformation by pure bending without stretching (see Novozhilov (Reference Novozhilov1959, pp. 25–26) for a discussion).

For scaling purposes, we approximate $A_1$ and $A_2$ as their values for an undeformed spherical surface of radius $R$, where (to recall) $R$ is the radius of the midsurface of the plate. Thus we have $A_1\approx R$ and $A_2\approx R\sin \theta$. Moreover, we neglect the terms $K_1 U$ in (3.7a,b), which our subsequent numerical solutions show to be small compared with ${A_1}^{-1}\partial _{\theta } W$. We then obtain

(3.8a,b)\begin{equation} \kappa_1\approx{-} \frac{1}{R^2}\partial^2_{\theta\theta} W, \quad \kappa_2\approx{-} \frac{\cot\theta}{R^2}\partial_{\theta} W. \end{equation}

Furthermore, (3.4) takes the form

(3.9)\begin{equation} F_{int} = \frac{1}{R} (M_1 - M_2)\cot\theta + \frac{1}{R}\partial_{\theta} M_1. \end{equation}

Now, using the scaling $(1/R) \partial _{\theta }\sim 1/l_b$, we find

(3.10)\begin{equation} M_1, M_2, M_1 - M_2\sim \frac{\eta_1 h^3V}{l_b^2} \left\langle1, \frac{l_b \cot\theta}{R}\right\rangle, \end{equation}

where the quantities inside $\langle \rangle$ indicate two terms with different scalings. Ignoring the small difference between $R$ and $R_0$ and choosing $\theta _t$ as a representative value of $\theta$, we define a dimensionless ‘sphericity number’

(3.11)\begin{equation} \varSigma = \frac{l_b}{R_0}\cot\theta_t. \end{equation}

Equations (3.9) and (3.10) now imply

(3.12)\begin{equation} F_{int} \sim \eta_0 V~f_1(St ,\varSigma) , \end{equation}

where

(3.13)\begin{equation} St = \frac{\eta_1}{\eta_0} \left( \frac{h}{l_b}\right)^3 \end{equation}

is a dimensionless ‘flexural stiffness’ and $f_1$ is an unknown function. Equation (3.13) is identical to the flexural stiffness for subduction of initially flat plates in two- and three-dimensional Cartesian geometries (Ribe Reference Ribe2010; Li & Ribe Reference Li and Ribe2012).

Having obtained the scalings for the three forces in (3.2), (3.3) and (3.12), we proceed to determine the scaling of the sinking speed $V$. Balancing $F_b$ and $F_{ext}$ yields the Stokes velocity scale

(3.14)\begin{equation} V_{Stokes} = \frac{l h g\delta\rho}{\eta_0}. \end{equation}

We expect the normalised sinking speed $V/V_{Stokes}$ to be a function of the ratio $F_{int}/F_{ext}$ and of the two angles $\theta _t$ (plate radius) and $\varphi _s$ (slab dip) whose purely geometrical influence cannot be captured by scaling analysis. We therefore obtain a scaling law of the general form

(3.15)\begin{equation} \frac{V}{V_{Stokes}}= f_2(St, \varSigma, \theta_t, \varphi_s) \equiv f_2\left( \frac{\eta_1}{\eta_0} \frac{h^3}{l_b^3}, \frac{l_b}{R_0}\cot\theta_t, \theta_t, \varphi_s \right) , \end{equation}

where $f_2$ is an unknown function. The ratio $d/h$ does not appear on the right-hand side of (3.15) because its only effect is to modify the bending length slightly (Ribe Reference Ribe2010). By introducing the intermediate variable $l_b$, our scaling analysis has succeeded in reducing the function (3.1) of six dimensionless arguments to (3.15), which involves only four arguments.

In the sequel we shall have occasion to compare our numerical predictions with those for a reference ‘flat-Earth’ limit in which spherical effects are absent. This limit corresponds to $\varSigma \to 0$, and can be achieved in two ways. The first is to note that $\varSigma$ may be written as

(3.16)\begin{equation} \varSigma = \frac{l_b}{L}\epsilon\cot\epsilon \quad \mathrm{where} \ \epsilon = \frac{L}{R_0} \end{equation}

and $L$ is the arcwise radius of the trench. Because $\lim _{\epsilon \to 0}\epsilon \cot \epsilon = 1$, we obtain $\varSigma \to 0$ in the double limit $L/R_0\to 0$, $l_b/L\to 0$. However, this limit is not easily accessible numerically, and so we define the flat-Earth limit in a different way as $\theta _t = {\rm \pi}/2$. This definition may not be immediately obvious, because it corresponds to a plate in the form of a complete hemisphere. However, the analogy becomes clearer when we consider that for a hemispherical plate the principal normal vector to the trench is perpendicular to the shell. The two-dimensional Cartesian (‘true flat-Earth’) case (Ribe Reference Ribe2010) is then recovered smoothly in the limit where the shell is vanishingly thin compared with the radius of curvature of the trench, which is approximately the radius of the Earth. The situation is fundamentally different for a spherical cap with $\theta _t < {\rm \pi}/2$, for which the principal normal has a component parallel to the shell's midsurface. We show in § 5.2 and figure 4 that the sinking speed $V(St)$ in the hemispherical limit $\theta _t = {\rm \pi}/2$ agrees closely with the two-dimensional Cartesian prediction for values of the flexural stiffness $St$ relevant to the Earth. We may therefore use $\theta _t={\rm \pi} /2$ as a physically sound proxy for the flat-Earth limit. The function (3.15) then simplifies to

(3.17)\begin{equation} \left(\frac{V}{V_{Stokes}} \right)_{flat} = f_2\left (St, 0, \frac{\rm \pi}{2}, \varphi_s\right )\equiv f_3(St,\varphi_s). \end{equation}

4. Boundary-element method

Our starting point is the general boundary-integral representation for the flow generated by a buoyant drop with density excess $\delta \rho$ and viscosity $\eta _1$ immersed in another fluid with viscosity $\eta _0$. Let $V_1$ denote the drop, $V_0$ denote the external fluid and $S$ denote the interface between them. Moreover, let $u_j^{(m)}(\boldsymbol {x})$ be the velocity in volume $V_m$. Then the flow inside and outside the drop is governed by the integral equation (Pozrikidis Reference Pozrikidis1990; Manga & Stone Reference Manga and Stone1993)

(4.1)\begin{align} \chi_0(\boldsymbol{x}_0) u_j^{(0)}(\boldsymbol{x}_0) + \gamma \chi_1(\boldsymbol{x}_0) u_j^{(1)}(\boldsymbol{x}_0) &={-}\frac{\delta\rho}{\eta_0}\int_S (\boldsymbol{g}(\boldsymbol{x})\boldsymbol{\cdot} \boldsymbol{x}){\mathsf{G}}_{ij}(\boldsymbol{x}, \boldsymbol{x}_0) n_i(\boldsymbol{x})\,\mathrm{d}S(\boldsymbol{x})\nonumber\\ &\quad + (1-\gamma)\int_S u_i(\boldsymbol{x}) {\mathsf{T}}_{ijk}(\boldsymbol{x}, \boldsymbol{x}_0) n_k(\boldsymbol{x})\,\mathrm{d}S(\boldsymbol{x}). \end{align}

Here, ${\mathsf{G}}_{ij}$ and ${\mathsf{T}}_{ijk}$ are Green functions for the velocity and stress, respectively, at the point $\boldsymbol {x}_0$ generated by a point force acting at $\boldsymbol {x}$. Also, $\chi _0(\boldsymbol {x}_0) = 1$ if $\boldsymbol {x}_0$ is in $V_0$, $1/2$ if $\boldsymbol {x}_0$ is right on $S$, and 0 if $\boldsymbol {x}_0$ is in $V_1$. $\chi _1(\boldsymbol {x}_0)$ is defined similarly but with the subscripts 0 and 1 interchanged. The velocity $\boldsymbol {u} (\boldsymbol {x}_0)$ on $S$ satisfies an integral equation obtained from (4.1) by setting $\chi _0 = \chi _1 = 1/2$ and applying the velocity matching condition $u_j^{(0)}(\boldsymbol {x}_0) = u_j^{(1)}(\boldsymbol {x}_0) = u_j(\boldsymbol {x}_0)$, yielding

(4.2)\begin{align} \frac{1}{2} (1 + \gamma) u_j(\boldsymbol{x}_0) &={-}\frac{\delta\rho}{\eta_0}\int_S (\boldsymbol{g}\boldsymbol{\cdot}\boldsymbol{x}){\mathsf{G}}_{ij}(\boldsymbol{x}, \boldsymbol{x}_0) n_i(\boldsymbol{x})\,\mathrm{d}S(\boldsymbol{x})\nonumber\\ &\quad + (1-\gamma)\int_S u_i(\boldsymbol{x}) {\mathsf{T}}_{ijk}(\boldsymbol{x}, \boldsymbol{x}_0) n_k(\boldsymbol{x})\,\mathrm{d}S(\boldsymbol{x}). \end{align}

Next, we subtract the singularities of the two integrands in (4.2) following the procedure outlined in § 6.4 of Pozrikidis (Reference Pozrikidis1992), which yields

(4.3a)\begin{equation} \int_S (\boldsymbol{g}\boldsymbol{\cdot}\boldsymbol{x}){\mathsf{G}}_{ij}(\boldsymbol{x}, \boldsymbol{x}_0) n_i(\boldsymbol{x})\,\mathrm{d} S(\boldsymbol{x}) = \int_S [\boldsymbol{g}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{x} - \boldsymbol{g}(\boldsymbol{x}_0) \boldsymbol{\cdot} \boldsymbol{x}_0]{\mathsf{G}}_{ij}(\boldsymbol{x}, \boldsymbol{x}_0) n_i(\boldsymbol{x})\,\mathrm{d} S(\boldsymbol{x}), \end{equation}

and

(4.3b)\begin{align} &\int_S u_i(\boldsymbol{x}){\mathsf{T}}_{ijk}(\boldsymbol{x}, \boldsymbol{x}_0) n_k(\boldsymbol{x})\,\mathrm{d} S(\boldsymbol{x})\nonumber\\ &\quad={-} \tfrac{1}{2} u_j (\boldsymbol{x}_0) + \int_S [u_i(\boldsymbol{x}) - u_i(\boldsymbol{x}_0)] {\mathsf{T}}_{ijk}(\boldsymbol{x}, \boldsymbol{x}_0) n_k(\boldsymbol{x})\,\mathrm{d} S(\boldsymbol{x}). \end{align}

Substituting (4.3b) into (4.2) we obtain

(4.4)\begin{align} u_j(\boldsymbol{x}_0) &={-} \frac{\delta\rho}{\eta_0}\int_S [\boldsymbol{g}(\boldsymbol{x}) \boldsymbol{\cdot}\boldsymbol{x} - \boldsymbol{g}(\boldsymbol{x}_0) \boldsymbol{\cdot} \boldsymbol{x}_0 ] {\mathsf{G}}_{ij}(\boldsymbol{x}, \boldsymbol{x}_0) n_i (\boldsymbol{x}) \,\mathrm{d}S(\boldsymbol{x})\nonumber\\ &\quad + (1 -\gamma)\int_S [u_i(\boldsymbol{x}) - u_i(\boldsymbol{x}_0)] {\mathsf{T}}_{ijk} (\boldsymbol{x}, \boldsymbol{x}_0) n_k(\boldsymbol{x})\,\mathrm{d} S(\boldsymbol{x}). \end{align}

We now set $\boldsymbol {g} = - g\boldsymbol {e}_r$ and non-dimensionalise all lengths by $R_0$ and all velocities by $g\delta \rho R_0^2/\eta _0$ (we do not non-dimensionalise using the velocity scale $V_{Stokes}$ because it contains the variable slab length $l$). Equation (4.4) then takes the form

(4.5a)\begin{gather} u_j(\boldsymbol{x}_0) ={-} \mathcal{S}_j(\boldsymbol{x}_0)+ (1- \gamma) \mathcal{D}_j(\boldsymbol{x}_0) , \end{gather}

(4.5b)\begin{gather}\mathcal{S}_j(\boldsymbol{x}_0) ={-} \int_S [\boldsymbol{e}_r(\boldsymbol{x}) \boldsymbol{\cdot} \boldsymbol{x} - \boldsymbol{e}_r(\boldsymbol{x}_0) \boldsymbol{\cdot}\boldsymbol{x}_0] n_i(\boldsymbol{x}) {\mathsf{G}}_{ij}(\boldsymbol{x}, \boldsymbol{x}_0)\,\mathrm{d}S(\boldsymbol{x}), \end{gather}

(4.5c)\begin{gather}\mathcal{D}_j(\boldsymbol{x}_0) = \int_S [u_i(\boldsymbol{x}) - u_i(\boldsymbol{x}_0)] {\mathsf{T}}_{ijk} (\boldsymbol{x}, \boldsymbol{x}_0) n_k(\boldsymbol{x})\,\mathrm{d} S(\boldsymbol{x}), \end{gather}

where all variables are now dimensionless. The quantity $\mathcal {S}_j(\boldsymbol {x}_0)$ is called the single-layer potential, and $\mathcal {D}_j(\boldsymbol {x}_0)$ is the double-layer potential.

The Green functions ${\mathsf{G}}_{ij}$ and ${\mathsf{T}}_{ijk}$ for flow inside a sphere are most readily derived by splitting the point force into a radial component normal to the sphere's surface and a transverse component tangential to it (Appendix A). It then proves convenient to work with mixed-index Green functions $\hat {{\mathsf{G}}}_{i\alpha }(\boldsymbol {x}, \boldsymbol {x}_0)$ and $\hat {{\mathsf{T}}}_{i\alpha k}(\boldsymbol {x}, \boldsymbol {x}_0)$, where the index $\alpha$ can take on either of the two values $r$ or $\theta$ according to whether the unit force points in the radial or the transverse direction at the location $\boldsymbol {x}_0$. Quantities involving mixed indices are indicated by a superposed hat. The key advantage of replacing the Cartesian index $j$ by the spherical polar index $\alpha$ is that, depending on the latter's value, we need use only the Green function for either a radial force or a transverse one to calculate $\hat {\boldsymbol{\mathsf{G}}}$ and hence the corresponding component of the single- and double-layer potentials.

Consider first the single-layer potential (4.5b), which we rewrite as

(4.6)\begin{equation} \hat{\mathcal{S}}_\alpha(\boldsymbol{x}_0)=\int_D (r_0-r)\hat{{\mathsf{G}}}_{i\alpha}(\boldsymbol{x},\boldsymbol{x}_0)n_i(\boldsymbol{x})\, \text{d}S(\boldsymbol{x}). \end{equation}

Next, we specialise (4.6) to axisymmetric flow by integrating over the azimuthal angle $\phi$. Following the notation of Pozrikidis (Reference Pozrikidis1992), we introduce cylindrical polar coordinates $(x,\sigma,\phi )$ such that the cylindrical axis $\sigma = 0$ corresponds to the Cartesian $x$-axis. We then set $\mathrm {d} S = \sigma \,\mathrm {d}\phi \,\mathrm {d} l$, where $\mathrm {d} l$ is the differential arclength of the trace of the contour $C$ in any azimuthal plane. We also note that the Cartesian components of the normal vector $\boldsymbol {n}$ are $(n_x, n_{\sigma }\cos \phi, n_{\sigma }\sin \phi )$. To perform the integration, we suppose for definiteness that the field point $\boldsymbol {x}_0$ lies in the $x$–$y$ plane where $\phi _0 = 0$. The single-layer potential can then be written as

(4.7)\begin{equation} \hat{\mathcal{S}}_\alpha(\boldsymbol{x}_0)=\int_C (r_0-r)\hat{\mathcal{M}}_{\alpha\beta}(\boldsymbol{x},\boldsymbol{x}_0)n_\beta(\boldsymbol{x})\, \text{d}l(\boldsymbol{x}), \end{equation}

where

(4.8) \begin{equation} \begin{pmatrix} \hat{\mathcal{M}}_{rx} & \hat{\mathcal{M}}_{r\sigma} \\ \hat{\mathcal{M}}_{\theta x} & \hat{\mathcal{M}}_{\theta\sigma} \end{pmatrix}= \sigma\int_0^{2{\rm \pi}}\begin{pmatrix} \hat{{\mathsf{G}}}_{xr} & \hat{{\mathsf{G}}}_{yr}\cos\phi+\hat{{\mathsf{G}}}_{zr}\sin\phi \\ \hat{{\mathsf{G}}}_{x\theta} & \hat{{\mathsf{G}}}_{y\theta}\cos\phi+\hat{{\mathsf{G}}}_{z\theta}\sin\phi \end{pmatrix}\text{d}\phi. \end{equation}

Note that only the radial Green function is necessary to calculate $\hat {\mathcal {M}}_{rx}$ and $\hat {\mathcal {M}}_{r\sigma }$, and hence the radial component $\hat {\mathcal {S}}_r$ of the single-layer potential. The situation for the transverse component is similar.

The procedure for the double-layer potential is analogous. The result is

(4.9)\begin{equation} \hat{\mathcal{D}}_\alpha(\boldsymbol{x}_0) =\int_C \left[ u_\beta(\boldsymbol{x})\hat{\mathcal{Q}}_{\alpha\beta \gamma}(\boldsymbol{x},\boldsymbol{x}_0) - u_\beta(\boldsymbol{x}_0)\hat{\mathcal{P}}_{\alpha\beta \gamma}(\boldsymbol{x},\boldsymbol{x}_0)\right] n_\gamma(\boldsymbol{x})\, \text{d}l(\boldsymbol{x}), \end{equation}

where the indices $\beta$ and $\gamma$ take cylindrical values $\{x,\sigma \}$ and the tensors $\hat {\mathcal {P}}$ and $\hat {\mathcal {Q}}$ are defined as

(4.10a)\begin{gather} \begin{pmatrix} \hat{\mathcal{P}}_{\alpha xx} & \hat{\mathcal{P}}_{\alpha x\sigma} \\ \hat{\mathcal{P}}_{\alpha \sigma x} & \hat{\mathcal{P}}_{\alpha \sigma\sigma} \end{pmatrix}= \sigma\int_0^{2{\rm \pi}}\begin{pmatrix} \hat{{\mathsf{T}}}_{x\alpha x} & \hat{{\mathsf{T}}}_{y\alpha x}\cos\phi+\hat{{\mathsf{T}}}_{z\alpha x}\sin\phi \\ \hat{{\mathsf{T}}}_{x\alpha y} & \hat{{\mathsf{T}}}_{y\alpha y}\cos\phi+\hat{{\mathsf{T}}}_{z\alpha y}\sin\phi \end{pmatrix}\text{d}\phi, \end{gather}

(4.10b)\begin{gather}\hat{\mathcal{Q}}_{\alpha xx} = \hat{\mathcal{P}}_{\alpha x x}, \quad \hat{\mathcal{Q}}_{\alpha x\sigma}= \hat{\mathcal{Q}}_{\alpha \sigma x} = \hat{\mathcal{P}}_{\alpha x\sigma}, \end{gather}

(4.10c)\begin{gather}\hat{\mathcal{Q}}_{\alpha \sigma\sigma} = \sigma\int_0^{2{\rm \pi}} \left( \hat{{\mathsf{T}}}_{y\alpha y}\cos^2\phi +\hat{{\mathsf{T}}}_{z\alpha z}\sin^2\phi+\hat{{\mathsf{T}}}_{y\alpha z}\sin2\phi\right)\mathrm{d}\phi. \end{gather}

Again, only the radial (transverse) Green function is needed to calculate the radial (transverse) component of the double-layer potential at $\boldsymbol {x}_0$.

To evaluate the single- and double-layer integrals (4.7) and (4.9), we discretise the contour $C$ using three-node curved elements $C_n$ $(n = 1, 2,\ldots, N)$, over each of which $\boldsymbol {x}$ and $\boldsymbol {u}$ vary as

(4.11a,b)\begin{equation} \boldsymbol{x}(\xi) = \sum_{m=1}^3\phi_m(\xi) \boldsymbol{x}_m, \quad \boldsymbol{u}(\xi) = \sum_{m=1}^3\phi_m(\xi) \boldsymbol{u}_m, \end{equation}

where $\boldsymbol {x}_m$ are the (known) nodal coordinates, $\boldsymbol {u}_m$ are the (unknown) nodal velocities and $\phi _m(\xi )$ are quadratic basis functions defined on a master element $\xi \in [-1,1]$. Substitution of (4.11a,b) into (4.7) and (4.9) with $\boldsymbol{x}_0\in C$ transforms the integrals over $C$ into sums of integrals over the elements $C_n$, each of which is mapped onto the interval $\xi \in [-1,1]$ and evaluated using 6-point Gauss–Legendre quadrature. The resulting system of $4N + 2$ coupled linear equations was solved using lower–upper (LU) decomposition and back substitution (Press et al. Reference Press, Teukolsky, Vetterling and Flannery1992), yielding the nodal velocities $\boldsymbol{u}_m$ with fourth-order accuracy. The correctness of the instantaneous solution for the velocity was verified against an analytical solution for two concentric spheres with different viscosities (Appendix C). For the few time-dependent cases we considered, the positions of all material points $\boldsymbol {x}\in C$ were advanced at each time step using explicit Euler stepping.

5. Instantaneous solutions

5.1. General features

Because inertia is negligible in planetary mantles, the evolutionary history of a subducting shell is nothing more than a sequence of quasi-static configurations whose dynamics is determined entirely by the shell's instantaneous shape. It therefore makes sense first to study the quasi-static dynamics of a shell as a function of its viscosity ratio and shape, without the added complexity of the (purely kinematic) time evolution. Thus, instead of regarding the model parameters shown in figure 1 ($\theta _t$, $\theta _s$, $\varphi _s$ etc.) as initial values for a time-dependent simulation, we treat them as free geometrical parameters that can be varied to represent a wide range of different shell shapes at some arbitrary instant in time.

To begin, we examine the instantaneous flow of a reference shell with the parameters $\theta _t = 30^{\circ }$, $\theta _s = 36^{\circ }$, $\varphi _s = 45^{\circ }$, $h/R_0 = 0.0157$, $d/h = 0.3$ and $\gamma = 100$. Figure 2(a) shows the velocity on the shell's midsurface, calculated by averaging the velocities on the upper and lower surfaces of the shell adjoining each midsurface point. The slab moves downward and backward, while the nearly horizontal velocity within the plate is of much smaller magnitude. Figure 2(b) shows the radial velocity as a function of arclength along the midsurface. The most interesting feature is an interval of upward radial velocity $s/R_0\in [0.49, 0.53]$. This corresponds to the so-called ‘flexural bulge’ of the ocean floor that is observed seaward of many trenches on Earth.

Figure 2.

Dynamics of an illustrative instantaneous subduction solution with $\theta _t = 30^{\circ }$, $\theta _s = 36^{\circ }$, $h/R_0 = 0.0157$, $d/h = 0.3$, $\varphi _s = 45^{\circ }$ and $\gamma = 100$. (a) Geometry of the shell with midsurface velocity vectors. The longest arrow has magnitude $0.299 h^2 g\delta \rho /\eta _0$. (b) Radial velocity $u_r$. The label ‘fb’ indicates a region of upward radial velocity corresponding to flexural bulging seaward of the trench. (c) Bending rate $\dot {K}$ and (d) stretching rate $\dot {E}$ as functions of arclength $s$ along the midsurface. (e) Rates of viscous dissipation of energy due to bending (red) and stretching (blue) as functions of $s$.

To proceed further with our analysis, we now exploit the fact that the shell is a thin object whose rate of deformation can be described as a combination of bending and stretching. According to the thin-shell constitutive relation (3.6a,b), the rate of bending of a line of constant longitude on the shell's midsurface is measured by the effective bending rate

(5.1)\begin{equation} \dot{K} ={-} \left(\kappa_1 + \tfrac{1}{2}\kappa_2\right), \end{equation}

where $\kappa _1$ ad $\kappa _2$ are defined by (3.7a,b). Similarly, the effective stretching rate of the midsurface is

(5.2)\begin{equation} \dot{E} = \epsilon_1 + \tfrac{1}{2}\epsilon_2, \end{equation}

where

(5.3a,b)\begin{equation} \epsilon_1 = \frac{1}{A_1}\partial_{\theta} U - K_1 W,\quad \epsilon_2 = \frac{1}{A_1 A_2} U\partial_{\theta} A_2 - K_2 W \end{equation}

and $K_2$ is the second principal curvature of the midsurface. Figures 2(c) and 2(d) show $\dot K(s)$ and $\dot E(s)$ for the reference shell. Significant bending is confined to a boundary layer of arcwise extent $\approx 0.14 R_0$ adjoining the edge of the shell. The inner (plateward) half of the boundary layer has $\dot {K} < 0$, which corresponds to counterclockwise bending. The outer half experiences clockwise bending ($\dot {K} > 0$) due to the stresses applied by the outer fluid flowing clockwise around the free edge of the slab. Turning to stretching, we see that both stretching ($\dot {E} > 0$) and compression ($\dot {E} < 0$) are significant in the boundary layer. Finally, the plate outside the boundary layer is under nearly uniform compression.

Because the bending rate $\dot {K}$ and the stretching rate $\dot {E}$ have different units, the best way to compare them is in terms of the rates of viscous dissipation of energy $\varPhi$ associated with deformation by bending ($\varPhi _b$) and stretching ($\varPhi _s)$. For axisymmetric flow, the explicit expressions for these dissipation rates per unit area of the midsurface are (Novozhilov Reference Novozhilov1959, (9.12))

(5.4a,b)\begin{equation} \varPhi_b = \tfrac{1}{6}\eta_1 h^3\left[ (\kappa_1 + \kappa_2)^2 - \kappa_1 \kappa_2\right], \quad \varPhi_s = 2\eta_1 h \left[(\epsilon_1 + \epsilon_2)^2 - \epsilon_1\epsilon_2\right]. \end{equation}

Figure 2(e) shows $\varPhi _b(s)$ (red line) and $\varPhi _s(s)$ (blue line) for the reference shell. The inner half of the boundary layer is strongly dominated by bending, whereas the outer half shows a rough equipartition of dissipation between bending and stretching.

5.2. Bending boundary layer

The next step is to define more precisely the width of the bending boundary layer, which is just the bending length $l_b$ that we introduced in our scaling analysis (§ 3.2). Figure 3 shows the bending rate $\dot {K}(s)$ for the reference shell geometry with viscosity ratios $\gamma = 10^2$ (red), $10^3$ (green) and $10^4$ (blue). As shown for the case $\gamma = 10^2$, we define $l_b$ as the distance from the free end of the slab to the first zero of $\dot {K}$ plateward (leftward in the figure) from the position where $\dot {K}$ is a minimum. All else being equal, $l_b$ increases as the viscosity ratio $\gamma$ increases: a stiffer shell bends over a greater distance.

Figure 3.

Bending rate $\dot {K}$ as a function of arclength for a shell with $\theta _t = 30^{\circ }$, $\theta _s = 36^{\circ }$, $h/R_0 = 0.0157$, $d/h = 0.3$, $\varphi _s = 45^{\circ }$ and $\gamma = 10^2$ (red), $10^3$ (green) and $10^4$ (blue). The definition of the bending length $l_b$ is shown for $\gamma = 10^2$.

To provide a basis of comparison for our subsequent solutions, we now characterise in more detail the flat-Earth scaling law from (3.17). We performed BEM simulations with $h/R_0 = 0.0157$, $d/h = 0.3$, $\theta _t = 90^{\circ }$, $\varphi _s = 45^{\circ }$, $\theta _s - \theta _t\in [2^{\circ }, 4^{\circ }]$ and $\gamma \in [10^2, 10^{4.25}]$. The resulting values of $V/V_{Stokes}$ are shown as a function of the flexural stiffness in figure 4. The points collapse onto a universal curve with two limits. For $St\ll 1$, the slope of the curve approaches zero on logarithmic axes, indicating that the sinking speed is proportional to the Stokes speed (3.14) and therefore entirely controlled by the exterior viscosity $\eta _0$. The operative balance here is between the external force $F_{ext}$ and the buoyancy force $F_{b}$ given by (3.3) and (3.2) respectively. For $St > 10$, the slope is equal to $-1$, implying that the sinking speed

(5.5)\begin{equation} V\sim \frac{g\delta\rho h l}{\eta_0}\left(\frac{\eta_1 h^3}{\eta_0 l_b^3}\right)^{{-}1} \propto \frac{1}{\eta_1} \end{equation}

depends only on the shell's own viscosity $\eta _1$. This behaviour corresponds to a balance between the buoyancy force and the internal viscous force $F_{int}$.

Figure 4.

Normalised sinking speed $V/V_{{Stokes}}$ as a function of the flexural stiffness $St$ for a slab with dip $\varphi _s = 45^{\circ }$. Open circles: predictions of the axisymmetric BEM code in the flat-Earth limit $\theta _t = 90^{\circ }$, with $d/h = 0.3$ and various combinations of values of $\gamma$ and $l/h$. Solid line: polynomial fit to the axisymmetric BEM predictions. Filled circles: predictions of the two-dimensional Cartesian (‘true flat-Earth’) BEM code of Ribe (Reference Ribe2010).

To verify that the hemispherical-plate limit $\theta _t = {\rm \pi}/2$ does indeed correspond to flat-Earth behaviour, we show in figure 4 the normalised sinking speed $V(St)/V_{Stokes}$ predicted by the two-dimensional Cartesian BEM code of Ribe (Reference Ribe2010) (filled circles). Remarkably, the axisymmetric and Cartesian predictions of $V(St)/V_{Stokes}$ agree to within 10 %–15 % for $St\leq 3$. This is all the more surprising given that the plate in the axisymmetric case is motionless while the two-dimensional plate moves towards the trench. Because $St<1$ is the relevant range of flexural stiffnesses for terrestrial subduction zones (see table 1), figure 4 demonstrates that the case $\theta _t = {\rm \pi}/2$ is a suitable proxy for subduction on a flat Earth.

Table 1.

Dimensionless parameters for Pacific subduction zones.

5.3. Effect of sphericity

Turning now to solutions with sphericity ($\varSigma > 0$), we first calculate the ratio $V/V_{flat}$ of the sphericity-influenced sinking speed $V$ to the flat-Earth sinking speed $V_{flat}$. The results are shown in figure 5 as functions of $St$ and $\varSigma$ for three values of the arcwise diameter $D$ of the plate: (a) 2224 km ($\theta _t = 10^{\circ }$), (b) 3335 km ($\theta _t = 15^{\circ }$) and (c) 4447 km ($\theta _t = 20^{\circ }$). Each panel is based on 50 BEM solutions with five different values of $\theta _s - \theta _t$ and ten different values of $\gamma$. In almost all cases $V/V_{flat} < 1$, indicating that sphericity reduces the sinking speed of the slab by stiffening the shell. In general, sphericity has a greater influence for longer slab length, greater shell stiffness $St$, and smaller plate diameter $D$. The last of these effects implies that the influence of sphericity is strongest for plates in the form of shallow shells, i.e. shells whose deviation from a plane is small compared with their radii of curvature. For the values of $\theta _s - \theta _t$ and $\gamma$ used, the maximum influence of sphericity on the sinking speed (in figure 5a) is nearly a factor of 4.

Figure 5.

Slab sinking speed $V$ normalised to the flat-Earth limit (figure 4) as a function of the flexural stiffness $St$ and the sphericity number $\varSigma$, for three values of the arcwise diameter $D$ of the plate: (a) 2224 km, (b) 3335 km and (c) 4447 km. Each panel comprises 50 numerical solutions with five different values of $\theta _s - \theta _t$ (labelled in panel (a)) and ten viscosity ratios $\gamma \in [10^2, 10^{4.25}]$ increasing from left to right.

A feature of figure 5 worth noting is that the angular plate radius $\theta _t$ influences $V/V_{flat}$ in two ways: via its appearance in the definition of the sphericity number $\varSigma = (l_b/R_0)\cot \theta _t$, and via its absolute value. Comparison of panels (a)–(c) of figure 5 immediately shows that $V/V_{flat}$ cannot be described as a function of $St$ and $\varSigma$ alone. This means that sphericity necessarily introduces two dimensionless parameters ($\varSigma$ and $\theta _t$) beyond those required to describe flat-Earth subduction, as already anticipated in the scaling law (3.15).

Proceeding as we just did for the sinking speed, we now quantify the effect of sphericity on the longitudinal normal stress (hoop stress) $\sigma _{22}$. This quantity is important because compressional hoop stress is what causes longitudinal buckling of the slab. Although such buckling cannot occur in our axisymmetric model, it nevertheless makes sense to calculate the axisymmetric hoop stress as a basic state whose stability to longitudinal perturbations can be analysed (as we shall do in a subsequent study.)

In thin-shell theory (Novozhilov Reference Novozhilov1959), the fundamental quantity related to the hoop stress is the stress resultant

(5.6)\begin{equation} T_2 = \int^{h/2}_{{-}h/2}\sigma_{22} (1 - K_1 z)\,\mathrm{d} z = 4\eta_1 h\left(\epsilon_2 + \tfrac{1}{2}\epsilon_1\right) , \end{equation}

where $\epsilon _1$ and $\epsilon _2$ are defined by (5.3a,b). The notation $T_2$ is that of Novozhilov (Reference Novozhilov1959); the single subscript prevents confusion with the stress Green function $T_{ijk}$. Figure 6 shows $T_2$ as a function of arclength for the same model parameters as in figure 2. Except for a small interval near the trench, $T_2 < 0$ everywhere, indicating compressive hoop stress. The maximum absolute value of $T_2$ occurs at the leading end of the slab.

Figure 6.

Hoop stress resultant $T_2$ as a function of arclength for the same model parameters as in figure 2.

Figure 7 shows $T_2$, normalised to its flat-Earth value $T_{2flat}$, as a function of the flexural stiffness $St$ and the sphericity number $\varSigma$. Depending on the values of $St$ and $\varSigma$, $T_2/T_{2flat}$ can be either greater than or (not much) less than unity. The largest effect of sphericity in this case occurs for relatively low values of $St$ and $\varSigma$. This is opposite to the case of the normalised sinking speed $V/V_{flat}$ (figure 5), for which the largest effect of sphericity occurs for relatively high values of $St$ and $\varSigma$.

Figure 7.

Same as figure 5, but for the hoop stress resultant $T_2$ normalised to its flat-Earth value $T_{2flat}$.

6. Illustrative time-dependent solutions

Having examined the scaling of the instantaneous subduction dynamics for fixed shell geometries, we now consider the time evolution of the shape of the shell itself. As an illustration, figure 8 shows the evolving shape of a shell for viscosity ratios $\gamma = 10^2$ and $10^3$, starting from the initial condition shown in black. Numbers near the end of each slab are dimensionless times $\hat {t} = t h^2 g\delta \rho /R_0\eta _0$. In both panels, the subducting shell gradually ‘peels away’ from the free surface such that the trench undergoes retrograde (plateward) motion. Geophysicists call this ‘trench rollback’. Comparing figures 8(a) and 8(b), we see that subduction evolves more rapidly for the lower viscosity ratio, simply because a less viscous shell can bend more easily. A further difference between the two cases concerns the shape of the slab. For $\gamma = 10^2$, the curvature of the lower part of the slab at $\hat {t} = 0.66$ is opposite in sign to that of the upper part. This is because the slab is weak enough to be deformed by the mantle ‘wind’ flowing counterclockwise around its end. By contrast, the curvature of the slabs with $\gamma = 10^3$ is of the same sign everywhere because the slab is stiff enough to resist deformation by the mantle wind. In closing, we re-emphasise that these purely axisymmetric solutions will no doubt be unstable to small longitudinal perturbations.

Figure 8.

Time-dependent subduction of shells with viscosity ratios (a) $\gamma = 10^2$ and (b) $10^3$. The initial condition for both cases (black) has $\theta _t = 20^{\circ }$, $\theta _s = 22^{\circ }$, $\varphi _s = 30^{\circ }$, $h/R_0 = 0.0157$ and $d/h = 0.3$. Numbers near the end of each slab give the time in units of $R_0\eta _0/h^2 g\delta \rho$.

7. Geophysical application

To apply our results to subduction on Earth, we now estimate the influence of sphericity on the dynamics of selected terrestrial subduction zones. Figure 9 shows the locations (green) of the six subduction zones we have chosen, all located in the Pacific ocean basin: Tonga, Marianas (both involving subduction of the Pacific plate), Ryukyu (Philippine Sea plate), Cascadia (Juan de Fuca plate), Central America (Cocos plate) and Chile (Nazca plate). Because our model is axisymmetric, we began by determining the radius $\theta _t$ of an equivalent spherical cap having the same area as the plate in question. These areas are listed in column 2 of table 1. Next, for each subduction zone we chose the model parameters $h$, $\ell$, $\theta _s$ and $\varphi _s$ by averaging values given by Lallemand, Heuret & Boutelier (Reference Lallemand, Heuret and Boutelier2005) for different transects perpendicular to the trench. We then ran instantaneous BEM simulations for the given slab geometry to determine the sinking speed $V$ and the hoop stress resultant $T_2$ for ten viscosity ratios $\gamma \in [10^2, 10^{4.25}]$. Finally, we ran ‘flat-Earth’ simulations for $\theta _t = 90^{\circ }$ and the same values of $\theta _s - \theta _t$ and the other parameters to obtain ten values of the flat-Earth sinking speed $V_{flat}$ and the hoop stress resultant $T_{2flat}$. Appendix D gives further details of these calculations.

Figure 9.

Map of the Pacific Ocean (Mercator projection) showing the boundaries of the major plates (red lines). Abbreviations of the plate names are shown in red: $\text {CO} = \text {Cocos}$, $\text {JF} = \text {Juan de Fuca}$, $\text {NZ} = \text {Nazca}$, $\text {PA} = \text {Pacific}$, $\text {PS} = \text {Philippine}$ Sea. The subduction zones that we consider are indicated by green lines, and their names are shown in green.

Figure 10 shows $V/V_{flat}$ as a function of the viscosity ratio $\gamma$ for three of our six subduction zones. Apart from one point in panel (a) with $\gamma = 100$, $V/V_{flat} < 1$, indicating that sphericity decreases $V$ relative to its flat-Earth value by stiffening the shell. In all cases, the effect of sphericity is stronger (smaller $V/V_{flat}$) for larger viscosity ratios. Finally, for a given viscosity ratio the sphericity effect is weak for the Pacific plate, stronger for the Nazca plate and strongest for the Cocos plate. For the geometry of the Cocos plate, the sphericity effect is nearly a factor of 4 for the highest viscosity ratio $\gamma = 10^{4.25}$.

Figure 10.

Effect of sphericity on the sinking speed for selected terrestrial subduction zones as a function of the assumed plate/mantle viscosity ratio $\gamma$. (a) Marianas (Pacific plate), (b) Chile (Nazca plate) and (c) Central America (Cocos plate). Each panel shows the numerically predicted slab sinking speed $V$ normalised to the flat-Earth limit as a function of the viscosity ratio. The red horizontal bars indicate the plausible range of $\gamma$ for the Earth.

To draw conclusions from figure 10 we need an estimate of the effective viscosity ratio $\gamma$ for subducting slabs on Earth. Several lines of evidence converge to suggest that $\gamma$ is of the order of a few hundred. Houseman & Gubbins (Reference Houseman and Gubbins1997) presented finite-element models of the Tonga slab and concluded that the modelled deformation matched the observed if $\gamma \approx 200$. Billen, Gurnis & Simons (Reference Billen, Gurnis and Simons2003) used regional finite-element models of the Tonga–Kermadec subduction zone to conclude that the magnitude of the observed strain rate requires $\eta _1 < 3\times 10^{23}$ Pa s or, equivalently, $\gamma < 300$. On the basis of analogue laboratory experiments, Funiciello et al. (Reference Funiciello, Faccenna, Heuret, Lallemand, Di Giuseppe and Becker2008) estimated that $\gamma \in [150,500]$ is required to explain the ratios of trench migration speeds to subducting plate speeds observed on Earth. An almost identical estimate $\gamma \in [140,510]$ was obtained by Ribe (Reference Ribe2010), who compared observed minimum radii of curvature of subducted slabs with the predictions of two-dimensional numerical models. Capitanio & Morra (Reference Capitanio and Morra2012) compared numerical model predictions with the observed inverse correlation between slab dip and radius of curvature to conclude that $\gamma \approx 200$.

The red horizontal bars in figure 10 show the range of $\gamma$ that brackets the aforementioned estimates. The maximum effect of sphericity occurs at the high end of the range ($\gamma = 510$), and is $1- V/V_{flat} = 7\,\%$ for Marianas, 20 % for Chile and 34 % for Central America.

After estimating $1- V/V_{flat}$ for our chosen subduction zones, we performed a similar calculation for the normalised hoop stress resultant $T_2/T_{2flat}$. The results of both calculations are summarised in table 1 for all six subduction zones and an assumed range of viscosity ratios $\gamma \in [140, 510]$. For the Tonga, Marianas and Chile subduction zones only the upper bound of $1- V/V_{flat}$ is given because calculation of the (small) lower bound is subject to numerical error. Also shown are the corresponding ranges of the flexural stiffness $St$ and the sphericity number $\varSigma$. While $St$ is small (${\leqslant }0.43$) for all six subduction zones, $\varSigma$ shows a strong inverse correlation with the plate size, as one would expect from the fact that $\varSigma \propto \cot \theta _t$. Finally, in all cases the effect of sphericity on the hoop stress resultant $T_2$ is much greater than the effect on the sinking speed $V$.