Alexander Chamolly

Neil M Ribe

Fluid mechanics of free subduction on a sphere, 1: The axisymmetric case

To understand how spherical geometry influences the dynamics of gravitydriven subduction of oceanic lithosphere on Earth, we study a simple model of a thin and dense axisymmetric shell of thickness h and viscosity η 1 sinking in a spherical body of fluid with radius R 0 and a lower viscosity η 0 . Using scaling analysis based on thin viscous shell theory, we identify a fundamental length scale, the 'bending length' l b , and two key dimensionless parameters that control the dynamics: the 'flexural stiffness' St = (η 1 /η 0 )(h/l b ) 3 and the 'sphericity number' Σ = (l b /R 0 ) cot θ t , where θ t is the angular radius of the subduction trench. To validate the scaling analysis, we obtain a suite of instantaneous numerical solutions using a boundary-element method based on new analytical point-force Green functions that satisfy free-slip boundary conditions on the sphere's surface. To isolate the effect of sphericity, we calculate the radial sinking speed V and the hoop stress resultant T 2 at the leading end of the subducted part of the shell, both normalised by their 'flat-Earth' values (i.e., for Σ = 0). For reasonable terrestrial values of η 1 /η 0 (≈ several hundred), sphericity has a modest effect on V , which is reduced by < 7% for large plates such as the Pacific plate and by up to 34% for smaller plates such as the Cocos and Philippine Sea plates. However, sphericity has a much greater effect on T 2 , increasing it by up to 64% for large plates and 240% for small plates. This result has important implications for the growth of longitudinal buckling instabilities in subducting spherical shells.

Introduction

Subduction of oceanic lithosphere is a major component of Earth's plate tectonic cycle: it is the main source of the buoyancy that drives mantle convection; it is the principal process responsible for recycling oceanic crust and volatile species like water back into the mantle; it is the main driver of long-term continental deformation; and it generates most of the great earthquakes and † Email address for correspondence: ribe@fast.u-psud.fr Abstract must not spill onto p.2 explosive volcanoes on Earth. Subduction occurs because oceanic lithosphere becomes denser as it cools moving away from the mid-ocean ridge where it formed. The lithosphere is therefore gravitationally unstable, and sinks into the mantle via a subcritical instability in which its negative buoyancy is sufficient to overcome its internal resistance to bending [START_REF] Mckenzie | The initiation of trenches: A finite amplitude instability[END_REF].

Since the classic 'corner flow' subduction model of [START_REF] Mckenzie | Speculations on the consequences and causes of plate motions[END_REF], hundreds of geodynamical models of subduction have been published, too numerous to cite here. The majority of these models have used two-or three-dimensional Cartesian geometry in which undeformed plates are flat. However, oceanic plates are doublycurved spherical shells, a fact that can be expected strongly to influence their mechanical behaviour. Whereas a flat plate supports a normal load by bending stresses alone, a doubly-curved shell can do so by a combination of bending stresses and in-plane 'membrane' stresses [START_REF] Audoly | Elasticity and Geometry: From Hair Curls to the Non-Linear Response of Shells[END_REF]. Moreover, curvature stiffens a shell and renders it more resistant to bending than a flat plate. This is because in order to to change the intrinsic curvature of a surface, additional energy needs to be expended on stretching. Familiar examples are a magazine or newspaper that one rolls up to swat a fly, or an orange peel that resists being flattened.

That terrestrial subduction occurs on a sphere is of course no secret in the geodynamics community, and the effect of sphericity has been the subject of numerous studies. Early models were purely geometrical, like the suggestion of [START_REF] Frank | Curvature of island arcs[END_REF] that the shape of the lithosphere in subduction zones resembles that of a circular dent surrounding a hole in a deformed ping-pong ball. [START_REF] Scholz | Buckling in island arcs[END_REF] and [START_REF] Bayly | Geometry of subducted plates and island arcs viewed as a buckling problem[END_REF] suggested that a subducting spherical shell should buckle along the strike of the trench to accommodate the reduction of the space available due to the spherical geometry as the shell penetrates deeper. [START_REF] Laravie | Geometry and lateral strain of subducted plates in island arcs[END_REF] proposed a geometric model that predicts different degrees of lateral strain in the subducted part of the shell depending on its dip and its radius of curvature at the surface. [START_REF] Schettino | Trench curvature and deformation of the subducting lithosphere[END_REF] reviewed the inadequacies of the pingpong-ball model and proposed an alternative kinematic model.

On the dynamical side, [START_REF] Tanimoto | Bending of a spherical lithosphere -axisymmetric case[END_REF][START_REF] Tanimoto | State of stress within a bending spherical shell and its implications for subducting lithosphere[END_REF] solved the equations for a normally loaded spherical elastic shell with negative buoyancy proportional to the (small) normal displacement, and concluded that the state of stress is strongly influenced by the spherical geometry. [START_REF] Mahadevan | Why subduction zones are curved[END_REF] used scaling analysis and numerical solutions for small-amplitude deformation of shallow spherical caps to investigate the causes of the curvature and segmentation of subduction zones. For the case of an elastic shell on a thicker elastic foundation, they found a scaling law for the wavelength of the edge instability ('dimpling') that occurs in response to a distributed radial body force. Finally, G. Morra and co-workers used the boundary-element method (BEM) to study large-amplitude subduction of viscous spherical shells, focusing on the curvature of island arcs [START_REF] Morra | Curvature of island arcs[END_REF], subduction of single plates in a mantle with depth-dependent viscosity, and interaction of multiple plates [START_REF] Morra | Spherical dynamic models of top-down tectonics[END_REF].

While the models discussed above focus narrowly on the subduction process itself, another approach is possible, namely time-dependent spherical thermal convection models in which subduction is but one aspect of the numerically predicted global circulation. In most such models, subduction occurs in an unrealistic twosided way, with two adjoining plates subducting together [START_REF] Coltice | What drives tectonic plates?[END_REF]). However, [START_REF] Schmeling | A benchmark comparison of spontaneous subduction models-Towards a free surface[END_REF] and [START_REF] Crameri | A free plate surface and weak oceanic crust produce single-sided subduction on Earth[END_REF] showed that realistic one-sided subduction could be obtained by adding to the top of the model domain a layer of low-viscosity fluid ('sticky air') to mimic a true free surface. From our point of view, the disadvantages of global circulation models are their high computational cost and relatively low spatial resolution, which typically prevent scaling laws from being sought.

A crucial question for our modelling concerns the effective long-term rheology of oceanic lithosphere, which is the key factor controlling its resistance to deformation. On short time scales typical of earthquakes and of the Earth's free oscillations, the lithosphere, like the rest of the solid Earth, behaves as an elastic medium that transmits shear and compressional waves. However, on long time scales characteristic of mantle convection (tens of Ma), the rheology of oceanic lithosphere is much more complicated [START_REF] Karato | Deformation of Earth Materials. An Introduction to the Rheology of Solid Earth[END_REF]. Experimental rock mechanics shows that the rheology of the lithosphere under conditions of strong bending is organised in three layers: brittle at shallow depths, elastic at depths around the neutral surface where the fibre stress vanishes, and ductile at deeper depths due to the higher temperatures there. Moreover, ductile behaviour can occur by three distinct mechanisms depending on temperature, pressure, stress and grain size: diffusion creep with a Newtonian rheology, dislocation creep with a generalised Newtonian (power-law) rheology, and low-temperature Peierls (exponential) creep. There is clearly much to be learned by including all these deformation mechanisms in a single realistic model (e.g. [START_REF] Bessat | Multiscale dynamics of the Tonga-Kermadec subduction zones[END_REF]). However, doing so tends to obscure the physical mechanisms at play and to prevent one from understanding scaling behaviour that a simpler model would reveal. Guided by these considerations, we have chosen to represent the rheology of the lithosphere by a constant viscosity that is much larger than that of the ambient mantle. In our view, this choice preserves the virtues of simplicity while embodying to lowest order the fact that the lithosphere is rheologically stiff relative to its surroundings.

In this paper, we present a simple model for the free (buoyancy-driven) subduction of an axisymmetric viscous spherical shell in an ambient fluid with a lower viscosity. Because the model involves two fluids separated by a sharp interface, it is amenable to solution by the boundary-element method (BEM). An important advantage of the BEM is that numerical solutions are highly accurate and quick to obtain, making possible the determination of clean quantitative scaling laws. Our model and results are novel in several ways. First, we focus on large-amplitude subduction, and take into account the self-consistent mechanical interaction between the highly deformed shell and the less viscous ambient mantle. Second, we identify a fundamental length scale, the 'bending length', that characterises the flexural response of the loaded shell. And third, we identify two key dimensionless parameters that govern the dynamics: a 'flexural stiffness' that determines which viscosity (shell or ambient) controls the sinking rate of the shell, and a 'sphericity number' that measures the importance of spherical geometry.

The paper is organised as follows. In § 2, we present the model problem and its geometrical and physical parameters. In § 3 we perform a dimensional analysis and a physical scaling analysis to reveal the length scale and dimensionless parameters that govern the dynamics. § 4 describes our implementation of the boundaryelement method. § 5 presents a suite of instantaneous BEM solutions, focusing on the effect of sphericity on the sinking speed of the shell and on the longitudinal normal stress (hoop stress) within it. § 6 presents illustrative time-dependent BEM solutions. Finally, § 7 estimates the magnitude of the sphericity effect for six subduction zones in the Pacific Ocean basin.

Model

The model envisions an axisymmetric viscous shell with density ρ 1 = ρ + δρ and viscosity η 1 immersed in a spherical viscous planet with density ρ 0 = ρ and viscosity η 0 (figure 1). The flow in both fluids is governed by the incompressible Stokes equations,

-∇p + η i ∇ 2 u + ρ i g = 0, ∇ • u = 0, (2.1)
where p is the pressure, u the velocity and ρ i 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 ≈ 0.54R 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 forces associated with its rotation, which is typically small compared to R 0 . 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 θ = θ t . The inner piece θ θ 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 θ t θ θ s is a downward-dipping ring of arcwise length l whose midsurface (dotted line in figure 1) is located at the radius

r(θ) = R(1 -bζ 3 -cζ 4 ), ζ = θ -θ t θ s -θ t (2.2)
where b and c are constants. The thickness of the ring is h everywhere outside the rounded rim just beyond θ = θ s . Following geophysical parlance, we shall call the cap θ < θ t the 'plate', the circle θ = θ t the 'trench' and the downward-dipping ring θ t θ θ s the 'slab' (these terms explain the choice of subscripts on θ t and θ s ). The mathematical form of Eq. (2.2) ensures that the local slope and local curvature (in the θ-direction) of the slab match those of the plate at the trench ζ = 0. The values of b and c are determined by imposing two additional constraints on Eq. (2.2). The first is that the dip ϕ of the leading end of the slab's midsurface relative to the local horizontal is a specified value ϕ 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 Focus on Fluids articles must not exceed this page length

are r (r 2 + (∂ θ r) 2 ) 1/2 ζ=1 = cos ϕ s , (2.3) - r 2 + 2(∂ θ r) 2 -r∂ 2 θθ r (r 2 + (∂ θ r) 2 ) 3/2 ζ=1 = - 1 R .
(2.4)

Eqs. (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, ϕ s and θ s -θ 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. [START_REF] Li | Dynamics of free subduction from 3-D boundary element modeling[END_REF] 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 [START_REF] Ribe | Bending mechanics and mode selection in free subduction: a thin-sheet analysis[END_REF], 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 neartrench 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.

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 geometry [START_REF] Ribe | Bending mechanics and mode selection in free subduction: a thin-sheet analysis[END_REF][START_REF] Li | Dynamics of free subduction from 3-D boundary element modeling[END_REF], 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 [START_REF] Ribe | Bending mechanics and mode selection in free subduction: a thin-sheet analysis[END_REF].

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 ϕ s , the angle θ t subtended by the plate, the mantle viscosity η 0 , the shell viscosity η 1 , and the buoyancy gδρ. The gravitational acceleration and the densities of the shell and the mantle appear only in the combination gδρ 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 Π-theorem [START_REF] Buckingham | On physically similar systems; illustrations of the use of dimensional equations[END_REF]) 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

V η 0 hlgδρ = f 0 h R 0 , h l , d h , η 1 η 0 , θ t , ϕ s (3.1)
where f 0 is an unknown function. The seven dimensionless groups could of course have been defined differently, but the choices in Eq. (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.

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

F b ∼ lhgδρ. (3.
2) 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 [START_REF] Ribe | Bending mechanics and mode selection in free subduction: a thin-sheet analysis[END_REF].

The viscous traction (force per unit area) applied to the shell by the outer fluid is σ ext ∼ η 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:

F ext ∼ η 0 V. (3.3)
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 [START_REF] Novozhilov | The Theory of Thin Shells[END_REF] together with the Stokes-Rayleigh analogy between incompressible elasticity and slow viscous flow [START_REF] Rayleigh | The Theory of Sound[END_REF]. According to this analogy, elastic shell theory can be transformed into its viscous equivalent by interpreting displacements as velocities and making the transformations E → 3η and σ → 1/2, where E is Young's modulus, η is the viscosity and σ is Poisson's ratio. From Eq. (7.8) of [START_REF] Novozhilov | The Theory of Thin Shells[END_REF], we have

F int ≡ N 1 = 1 A 1 A 2 [∂ θ (A 2 M 1 ) -M 2 ∂ θ A 2 ] , (3.4)
where

A 1 = r 2 + (∂ θ r) 2 1/2
and A 2 = r sin θ (3.5) are the Lamé parameters of the deformed axisymmetric shell and M 1 and M 2 are bending moments. The explicit expressions for these are

M 1 = ηh 3 3 κ 1 + 1 2 κ 2 , M 2 = ηh 3 3 κ 2 + 1 2 κ 1 (3.6)
where

κ 1 = - 1 A 1 ∂ θ 1 A 1 ∂ θ W + K 1 U , κ 2 = - 1 A 1 A 2 ∂ θ A 2 1 A 1 ∂ θ W + K 1 U (3.7)
and U and W are the velocity components tangential to and normal to the midsurface, respectively. The quantity K 1 in Eq. (3.7) is the curvature of a line of constant longitude; its explicit expression is the left-hand side of Eq. (2.4). Physically, the quantities -κ 1 and -κ 2 are the rates of change of curvature of the midsurface in the θ-and φ-directions, respectively, due to deformation by pure bending without stretching (see [START_REF] Novozhilov | The Theory of Thin Shells[END_REF] 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 ≈ R and A 2 ≈ R sin θ. Moreover, we neglect the terms K 1 U in Eq. (3.7), which our subsequent numerical solutions show to be small compared to A 1 -1 ∂ θ W . We then obtain

κ 1 ≈ - 1 R 2 ∂ 2 θθ W, κ 2 ≈ - cot θ R 2 ∂ θ W. (3.8)
Furthermore, Eq. (3.4) takes the form

F int = 1 R (M 1 -M 2 ) cot θ + 1 R ∂ θ M 1 . (3.9) Now using the scaling (1/R)∂ θ ∼ 1/l b , we find M 1 , M 2 , M 1 -M 2 ∼ η 1 h 3 V l 2 b 1, l b cot θ R , (3.10)
where the quantities inside indicate two terms with different scalings. Ignoring the small difference between R and R 0 and choosing θ t as a representative value of θ, we define a dimensionless 'sphericity number'

Σ = l b R 0 cot θ t .
(3.11)

Eqs. (3.9) and (3.10) now imply

F int ∼ η 0 V f 1 (St, Σ) (3.12)
where

St = η 1 η 0 h l b 3 (3.13)
is a dimensionless 'flexural stiffness' and f 1 is an unknown function. Eq. (3.13) is identical to the flexural stiffness for subduction of initially flat plates in two-and three-dimensional Cartesian geometry [START_REF] Ribe | Bending mechanics and mode selection in free subduction: a thin-sheet analysis[END_REF][START_REF] Li | Dynamics of free subduction from 3-D boundary element modeling[END_REF].

Having obtained the scalings for the three forces in Eqs. (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

V Stokes = lhgδρ η 0 . (3.14)
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 θ t (plate radius) and ϕ s (slab dip) whose purely geometrical influence cannot be captured by scaling analysis. We therefore obtain a scaling law of the general form

V V Stokes = f 2 (St, Σ, θ t , ϕ s ), ≡ f 2 η 1 η 0 h 3 l 3 b , l b R 0 cot θ t , θ t , ϕ s (3.15)
where f 2 is an unknown function. The ratio d/h does not appear on the right-hand side of Eq. (3.15) because its only effect is to modify the bending length slightly [START_REF] Ribe | Bending mechanics and mode selection in free subduction: a thin-sheet analysis[END_REF]. By introducing the intermediate variable l b , our scaling analysis has succeeded in reducing the function Eq. (3.1) of six dimensionless arguments to Eq. (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 Σ → 0, and can be achieved in two ways. The first is to note that Σ may be written as

Σ = l b L cot where = L R 0 (3.16)
and L is the arcwise radius of the trench. Because lim →0 cot = 1, we obtain Σ → 0 in the double limit L/R 0 → 0, l b /L → 0. However, this limit is not easily accessible numerically, and so we define the flat-Earth limit in a different way as θ t = π/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 [START_REF] Ribe | Bending mechanics and mode selection in free subduction: a thin-sheet analysis[END_REF]) is then recovered smoothly in the limit where the shell is vanishingly thin compared to 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 θ t < π/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 θ t = π/2 agrees closely with the two-dimensional Cartesian prediction for values of the flexural stiffness St relevant to the Earth. We may therefore use θ t = π/2 as a physically sound proxy for the flat-Earth limit. The function Eq. (3.15) then simplifies to

V V Stokes f lat = f 2 St, 0, π 2 , ϕ s ≡ f 3 (St, ϕ s ).
(3.17)

Boundary-element method

Our starting point is the general boundary-integral representation for the flow generated by a buoyant drop with density excess δρ and viscosity η 1 immersed in another fluid with viscosity η 0 . Let V 1 denote the drop, V 0 denote the external fluid, and S denote the interface between them. Moreover, let u (m) j (x) be the velocity in volume V m . Then the flow inside and outside the drop is governed by the integral equation [START_REF] Pozrikidis | The deformation of a liquid drop moving normal to a plane wall[END_REF][START_REF] Manga | Buoyancy-driven interaction between two deformable viscous drops[END_REF]

χ 0 (x 0 )u (0) j (x 0 ) + γχ 1 (x 0 )u (1) j (x 0 ) = - δρ η 0 S (g(x) • x)G ij (x, x 0 )n i (x)dS(x) + (1 -γ) S u i (x)T ijk (x, x 0 )n k (x)dS(x). (4.1)
Here G ij and T ijk are Green functions for the velocity and stress, respectively, at the point x 0 generated by a point force acting at x. Also, χ 0

(x 0 ) = 1 if x 0 is in V 0 , 1/2 if x 0 is right on S, and 0 if x 0 is in V 1 . χ 1 (x 0
) is defined similarly but with the subscripts 0 and 1 interchanged. The velocity u(x 0 ) on S satisfies an integral equation obtained from Eq. (4.1) by setting χ 0 = χ 1 = 1/2 and applying the velocity matching condition u (0)

j (x 0 ) = u (1) j (x 0 ) = u j (x 0 ), yielding 1 2 (1 + γ)u j (x 0 ) = - δρ η 0 S (g • x)G ij (x, x 0 )n i (x)dS(x) + (1 -γ) S u i (x)T ijk (x, x 0 )n k (x)dS(x). (4.2)
Next we subtract the singularities of the two integrands in Eq. (4.2) following the procedure outlined in §6.4 of [START_REF] Pozrikidis | Boundary Integral and Singularity Methods for Linearized Viscous Flow[END_REF], which yields

S (g • x)G ij (x, x 0 )n i (x)dS(x) = S [g(x) • x -g(x 0 ) • x 0 ]G ij (x, x 0 )n i (x)dS(x), (4.3a)
and

S u i (x)T ijk (x, x 0 )n k (x)dS(x) = - 1 2 u j (x 0 ) + S [u i (x) -u i (x 0 )]T ijk (x, x 0 )n k (x)dS(x). (4.3b)
Substituting Eq. (4.3b) into Eq. ( 4.2) we obtain

u j (x 0 ) = - δρ η 0 S [g(x) • x -g(x 0 ) • x 0 ]G ij (x, x 0 )n i (x)dS(x) + (1 -γ) S [u i (x) -u i (x 0 )]T ijk (x, x 0 )n k (x)dS(x). (4.4)
We now set g = -ge r and non-dimensionalise all lengths by R 0 and all velocities by gδρR 2 0 /η 0 (we do not non-dimensionalise using the velocity scale V Stokes because it contains the variable slab length l). Eq. (4.4) then takes the form

u j (x 0 ) = -S j (x 0 ) + (1 -γ)D j (x 0 ), (4.5a) S j (x 0 ) = - S [e r (x) • x -e r (x 0 ) • x 0 ]n i (x)G ij (x, x 0 )dS(x), (4.5b) D j (x 0 ) = S [u i (x) -u i (x 0 )]T ijk (x, x 0 )n k (x)dS(x), (4.5c)
where all variables are now dimensionless. The quantity S j (x) is called the singlelayer potential, and D j (x) is the double-layer potential.

The Green functions G ij and 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 Ĝiα (x, x 0 ) and Tiαk (x, x 0 ), where the index α can take on either of the two values r or θ according to whether the unit force points in the radial or the transverse direction at the location 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 α is that depending on the latter's value we need use only the Green function for either Rapids articles must not exceed this page length a radial force or a transverse one to calculate Ĝ and hence the corresponding component of the single-and double-layer potentials.

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

Ŝα (x 0 ) = D (r 0 -r) Ĝiα (x, x 0 )n i (x) dS(x). (4.6)
Next, we specialise Eq. (4.6) to axisymmetric flow by integrating over the azimuthal angle φ. Following the notation of [START_REF] Pozrikidis | Boundary Integral and Singularity Methods for Linearized Viscous Flow[END_REF], we introduce cylindrical polar coordinates (x, σ, φ) such that the cylindrical axis σ = 0 corresponds to the Cartesian x-axis. We then set dS = σdφdl, where dl 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 n are (n x , n σ cos φ, n σ sin φ). To perform the integration, we suppose for definiteness that the field point x 0 lies in the x-y plane where φ 0 = 0. The single-layer potential can then be written as To evaluate the single-and double-layer integrals Eq. (4.7) and Eq. (4.9), we discretise the contour C using three-node curved elements C n (n = 1, 2, ..., N ), over each of which x and u vary as (4.11) where x m are the (known) nodal coordinates, u m are the (unknown) nodal velocities, and φ m (ξ) are quadratic basis functions defined on a master element ξ ∈ [-1, 1]. Substitution of Eq. (4.11) into Eq. (4.7) and Eq. (4.9) with x 0 ∈ C transforms the integrals over C into sums of integrals over the elements C n , each of which is mapped onto the interval ξ ∈ [-1, 1] and evaluated using 6-point Gauss-Legendre quadrature. The resulting system of 4N + 2 coupled linear equations was solved using LU decomposition and back-substitution [START_REF] Press | Numerical Recipes in Fortran 77[END_REF], yielding the nodal velocities 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 x ∈ C were advanced at each time step using explicit Euler stepping.

Ŝα (x 0 ) = C (r 0 -r) Mαβ (x, x 0 )n β (x) dl(x), ( 4
x(ξ) = 3 m=1 φ m (ξ)x m , u(ξ) = 3 m=1 φ m (ξ)u m ,
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 are 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 (θ t , θ s , ϕ 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 θ t = 30 • , θ s = 36 • , ϕ s = 45 • , h/R 0 = 0.0157, d/h = 0.3, and γ = 100. Figure 2a 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 2b 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 ∈ [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.

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 Eq. (3.6), the rate of bending of a line of constant longitude on the shell's midsurface is measured by the effective bending rate

K = -κ 1 + 1 2 κ 2 (5.1)
where κ 1 ad κ 2 are defined by Eq. (3.7). Similarly, the effective stretching rate of the midsurface is

Ė = 1 + 1 2 2 (5.2)
where

1 = 1 A 1 ∂ θ U -K 1 W, 2 = 1 A 1 A 2 U ∂ θ A 2 -K 2 W (5.3)
and K 2 is the second principal curvature of the midsurface. Figures 2c and2d show K(s) and Ė(s) for the reference shell. Significant bending is confined to a

.

Φη /h g δρ stretching bending boundary layer of arcwise extent ≈ 0.14R 0 adjoining the edge of the shell. The inner (plateward) half of the boundary layer has K < 0, which corresponds to counterclockwise bending. The outer half experiences clockwise bending ( 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 ( Ė > 0) and compression ( Ė < 0) are significant in the boundary layer. Finally, the plate outside the boundary layer is under nearly uniform compression.

Because the bending rate K and the stretching rate Ė have different units, the best way to compare them is in terms of the rates of viscous dissipation of energy Φ associated with deformation by bending (Φ b ) and stretching (Φ s ). For axisymmetric flow, the explicit expressions for these dissipation rates per unit area of the midsurface are [START_REF] Novozhilov | The Theory of Thin Shells[END_REF], Eq. (9.12))

Φ b = 1 6 η 1 h 3 (κ 1 + κ 2 ) 2 -κ 1 κ 2 , Φ s = 2η 1 h ( 1 + 2 ) 2 -1 2 .
(5.4)

Figures 2e shows Φ b (s) (red line) and Φ 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.

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 K(s) for the reference shell geometry with viscosity ratios γ = 10 2 (red), 10 3 (green) and 10 4 (blue). As shown for the case γ = 10 2 , we define l b as the distance from the free end of the slab to the first zero of K plateward (leftward in the figure) from the position where K is a minimum. All else being equal, l b increases as the viscosity ratio γ increases: a stiffer shell bends over a greater distance.

To provide a basis of comparison for our subsequent solutions, we now characterise in more detail the flat-Earth scaling law from Eq. (3.17). We performed BEM simulations with h/R 0 = 0.0157, d/h = 0.3, 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 1, the slope of the curve approaches zero on logarithmic axes, indicating that the sinking speed is proportional to the Stokes speed Eq. (3.14) and therefore entirely controlled by the exterior viscosity η 0 . The operative balance here is between the external force F ext and the buoyancy force F b given by Eqs. (3.3) and (3.2) respectively. For St > 10, the slope is equal to -1, implying that the sinking speed

θ t = 90 • , ϕ s = 45 • , θ s -θ t ∈ [2 • , 4 • ] and γ ∈ [10 2 ,
V ∼ gδρhl η 0 η 1 h 3 η 0 l 3 b -1 ∝ 1 η 1 (5.5)
depends only on the shell's own viscosity η 1 . This behaviour corresponds to a balance between the buoyancy force and the internal viscous force F int .

To verify that the hemispherical-plate limit θ t = π/2 does indeed correspond to flat-Earth behaviour, we show in figure 4 (2010) (filled circles). Remarkably, the axisymmetric and Cartesian predictions of V (St)/V Stokes agree to within 10-15% for St 3. This is all the more surprising given that the plate in the axisymmetric case is motionless while the twodimensional 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 θ t = π/2 is a suitable proxy for subduction on a flat Earth.

Effect of sphericity

Turning now to solutions with sphericity (Σ > 0), we first calculate the ratio V /V f lat of the sphericity-influenced sinking speed V to the flat-Earth sinking speed V f lat . The results are shown in figure 5 In almost all cases V /V f lat < 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 to their radii of curvature. For the values of θ s -θ t and γ used, the maximum influence of sphericity on the sinking speed (in figure 5a) is nearly a factor of 4. A feature of figure 5 worth noting is that the angular plate radius θ t influences V /V f lat in two ways: via its appearance in the definition of the sphericity number Σ = (l b /R 0 ) cot θ t , and via its absolute value. Comparison of parts (a)-(c) of figure 5 immediately shows that V /V f lat cannot be described as a function of St and Σ alone. This means that sphericity necessarily introduces two dimensionless parameters (Σ and θ t ) beyond those required to describe flat-Earth subduction, as already anticipated in the scaling law Eq. (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) σ 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 [START_REF] Novozhilov | The Theory of Thin Shells[END_REF], the fundamental quantity related to the hoop stress is the stress resultant

T 2 = h/2 -h/2 σ 22 (1 -K 1 z)dz = 4η 1 h 2 + 1 2 1 (5.6)
where 1 and 2 are defined by Eq. ( 5.3). The notation T 2 is that of Novozhilov (1959); 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. and Σ, T 2 /T 2f lat 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 Σ. This is opposite to the case of the normalised sinking speed V /V f lat (figure 5), for which the largest effect of sphericity occurs for relatively high values of St and Σ.

Illustrative time-dependent solutions

Having examined the scaling of 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 γ = 10 2 and 10 3 , starting from the initial condition shown in black. Numbers near the end of each slab are dimensionless times t = th 2 gδρ/R 0 η 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 8a and 8b, 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 γ = 10 2 , the curvature of the lower part of the slab at 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 γ = 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.

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. model is axisymmetric, we began by determining the radius θ 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, , θ s and ϕ s by averaging values given by [START_REF] Lallemand | On the relationships between slab dip, backarc stress, upper plate absolute motion, and crustal nature in subduction zones[END_REF] 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 γ ∈ [10 2 , 10 4.25 ]. Finally, we ran 'flat-Earth' simulations for θ t = 90 • and the same values of θ s -θ t and the other parameters to obtain ten values of the flat-Earth sinking speed V f lat and the hoop stress resultant T 2f lat . Appendix D gives further details of these calculations.

Figure 10 shows V /V f lat as a function of the viscosity ratio γ for three of our six subduction zones. Apart from one point in part (a) with γ = 100,

V /V f lat < 1, Plate Area (km 2 ) Subduction zone St Σ 1 -V /V f lat T2/T 2f lat
PA 1.05 × 10 8 Tonga 0.13-0.37 0.12-0.13 0.065 1.33-1.37 PA 1.05 × 10 8 Marianas 0.12-0.33 0.12-0.13 0.069 1.45-1.64 NZ 1.61 × 10 7 Chile 0.062-0.17 0.48-0.52 0.20 1.29-1.44 PS 5.44 × 10 6 Ryukyu 0.16-0.43 0.56-0.62 0.11-0.33 1.91-2.39 CO 2.93 × 10 6 Central America 0.091-0.25 0.66-0.72 0.11-0.34 1.91-2.31 JF 2.56 × 10 5 Cascadia 0.029-0.083 2.3-2.5 0.12-0.33 1.69-2.19

Table 1: Dimensionless parameters for Pacific subduction zones. 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 f lat ) 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 γ = 10 4.25 .

To draw conclusions from figure 10 we need an estimate of the effective viscosity ratio γ for subducting slabs on Earth. Several lines of evidence converge to suggest that γ is on the order of a few hundred. [START_REF] Houseman | Deformation of subducted oceanic lithosphere[END_REF] presented finite-element models of the Tonga slab and concluded that the modelled deformation matched the observed if γ ≈ 200. [START_REF] Bessat | Multiscale dynamics of the Tonga-Kermadec subduction zones[END_REF] used regional finite-element models of the Tonga-Kermadec subduction zone to conclude that the magnitude of the observed strain rate requires η 1 < 3 × 10 23 Pa.s or, equivalently, γ < 300. On the basis of analog laboratory experiments, [START_REF] Funiciello | Trench migration, net rotation and slab-mantle coupling[END_REF] estimated that γ ∈ [150, 500] is required to explain the ratios of trench migration speeds to subducting plate speeds observed on Earth. An almost identical estimate γ ∈ [140, 510] was obtained by [START_REF] Ribe | Bending mechanics and mode selection in free subduction: a thin-sheet analysis[END_REF], who compared observed minimum radii of curvature of subducted slabs with the predictions of two-dimensional numerical models. [START_REF] Capitanio | The bending mechanics in a dynamic subduction system: Constraints from numerical modelling and global compilation analysis[END_REF] compared numerical model predictions with the observed inverse correlation between slab dip and radius of curvature to conclude that γ ≈ 200.

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

After estimating 1 -V /V f lat for our chosen subduction zones, we performed a similar calculation for the normalised hoop stress resultant T 2 /T 2f lat . The results of both calculations are summarised in table 1 for all six subduction zones and an assumed range of viscosity ratios γ ∈ [140,510]. For the Tonga, Marianas and Chile subduction zones only the upper bound of 1 -V /V f lat 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 Σ. While St is small ( 0.43) for all six subduction zones, Σ shows a strong inverse correlation with the plate size, as one would expect from the fact that Σ ∝ cot θ 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 .

Discussion

We begin by discussing further some of the limitations of our model, their motivations, and how they might be overcome. The first obvious limitation is the model's axisymmetry. One consequence of this is that the plate itself does not move, so that subduction occurs by trench rollback alone. This situation is not encountered on Earth, where subducting plates usually have (crudely speaking) a mid-ocean ridge on one side and a subduction zone on the other. The plate therefore can be pulled by the subducting slab ('slab pull' in geophysical parlance), causing it to move from the ridge towards the trench. Axisymmetry further implies that the trench is everywhere convex inward. While this is the case on the large scale, a glance at figure 9 shows that the subduction zones in the western Pacific ocean are composed of adjoining concave-inward segments. However, the purpose of the axisymmetry assumption is not geophysical realism, but rather to formulate the simplest possible model that allows us to quantify the influence of sphericity on key dynamical aspects of subduction. Once the essential goal of identifying the relevant length scales and dimensionless parameters has been attained, one can proceed to investigate three-dimensional subduction in more realistic geometries that embody the features noted above.

The second major simplification in our model is our neglect of an effectively inviscid core, which on Earth has a radius 0.54R 0 . Such a core acts as a freeslip boundary condition at the bottom of the mantle, and will have a significant quantitative influence on the flow generated by subducting slabs [START_REF] Morra | Earth curvature effects on subduction morphology: Modeling subduction in a spherical setting[END_REF][START_REF] Morra | Spherical dynamic models of top-down tectonics[END_REF]. We have two reasons for ignoring the core. First, its presence would introduce an additional length scale into the model, making it harder to determine the underlying scalings while providing no compensating physical insight. Second, because no Green functions for a spherical annulus bounded by two free-slip surfaces exist in the literature, a three-dimensional BEM approach would become unwieldy, requiring discretisation of the entire surfaces of the Earth and of the core. At this point one would probably be well advised to replace the BEM by more flexible approaches such as finite element or finite volume methods, which allow additional realistic aspects of the mantle (mineralogical phase transitions, radial variation of viscosity, etc.) to be incorporated easily. As a final remark, we emphasise that the quantities V /V f lat and T 2 /T 2f lat that we calculate in this study are normalised quantities with the flat-Earth limit in the denominator. It is likely that the presence of a core will influence the numerator and the denominator similarly. If so, then our estimates of the normalised quantities will remain valid even if a core is present.

Our results show that the effect of sphericity depends strongly on the size of the subducting plate in a way that defies naive intuition. At first glance, one might expect the effect of sphericity to be greatest for a large plate, whose midsurface differs more from a plane than that of a small or 'shallow' plate. However, the truth is exactly the opposite. Mathematically speaking, this is so because the dimensionless number that measures the effect of sphericity is Σ = (l b /R 0 ) cot θ t . For the same value of l b , therefore, Σ is greater for a small plate (small θ t ) than for a large one. Physically, one can understand this result by cutting two spherical shells from a child's plastic ball: a large one equal to a full hemisphere, and a small one in the form of a shallow spherical cap. Now balance each shell upside down on the point of an upright pencil, and deform the shell slightly by applying radially directed forces to opposite sides. You will then see that a given applied force produces a smaller deformation of the shallow cap than of the hemisphere. In other words, the hemisphere is 'floppy' while the shallow cap is 'stiff'. The same result applies to viscous shells if one replaces 'deformation' by 'rate of deformation'.

From a geophysical point of view, the main conclusion of this work is that sphericity affects the longitudinal normal stress resultant T 2 much more than the sinking speed V . This accords with the conclusion of [START_REF] Tanimoto | State of stress within a bending spherical shell and its implications for subducting lithosphere[END_REF] that the dominant effect of sphericity is on the state of stress in subducted lithosphere. We find that the effect of sphericity on V is rather modest, being at most about 7% for the largest (Pacific) plate and 34% for small plates (Philippine Sea, Cocos, Juan de Fuca) with θ t 12 • . By contrast, the effect on T 2 is up to 64% for the Pacific plate and up to 240% for the small plates. Our numerical solutions show that T 2 < 0 in the slab, which corresponds to compressional hoop stresses. These stresses arise because a slab subducting in a sphere must squeeze into an ever smaller space as it descends, due to the diminution with depth of the surface area of concentric spheres. While this diminution could in principle be accommodated by uniform thickening of the slab, it is more likely that the compressional hoop stresses drive longitudinal buckling instabilities with growth rates enhanced by sphericity [START_REF] Morra | Earth curvature effects on subduction morphology: Modeling subduction in a spherical setting[END_REF][START_REF] Morra | Spherical dynamic models of top-down tectonics[END_REF]. The implication is that the axisymmetric solutions we have discussed here will be unstable to small non-axisymmetric perturbations. A BEM investigation of three-dimensional spherical subduction with dynamic buckling is underway and will be reported separately.

Both authors thank E. Lauga for bringing us together. E. Calais kindly provided a file of Earth's tectonic plate boundaries. We thank G. Morra and an anonymous referee for constructive comments that helped us significantly to improve the manuscript. N.M.R. was supported by the Progamme National de Planétologie (PNP) of the Institut des Sciences de l'Univers (INSU) of the CNRS, co-funded by CNES.

(0, 0, 0)

(x, σ, ϕ) (x , σ , 0) f r -f/r θ (x , σ , 0) IM IM x -x x -x IM x σ
≡ r -1 0 cos θ0, σ IM 0 ≡ r -1 0 sin θ0, 0).
The force acting at the image point is -r -1 0 f . An arbitrary integration point has cylindrical coordinates (x, σ, φ) with the same values of (r0, θ0) as the field point.

A.2. Transverse point force

The Green functions for this case can be obtained by adapting the solution of [START_REF] Padmavathi | Motion inside a liquid sphere: internal singularities[END_REF] for flows inside an impermeable spherical boundary enclosing singularities and separating fluids having different viscosities. The limit relevant to our problem is when the outer viscosity is zero, corresponding to a free-slip surface.

We begin by writing the solution of the Stokes equations in terms of two scalar potentials {A, B}, as

u = ∇ × (∇ × (Ax)) + ∇ × (Bx) , (A 5a) p = p 0 + η ∂ ∂r r∇ 2 A , ( A 5b)
where r = √ x 2 + y 2 + z 2 and the potentials satisfy

∇ 4 A = 0, ∇ 2 B = 0. (A 6)
The condition of impermeability of the surface r = 1 is satisfied if

A = 0 on r = 1. (A 7)
The no-stress boundary condition, σ rθ = σ rφ = 0 is satisfied when

∂ 2 A ∂r 2 = 0, ∂ ∂r B r = 0 on r = 1. (A 8)
We therefore have a total of three boundary conditions that the two scalar fields A and B must satisfy.

The potentials for a transverse Stokeslet with strength f θ e θ,0 are

A 0 (x; x 0 ) = f θ 8πη x • e r,0 -r 0 r 0 (R -r) + r -r 0 x • e r,0 r x • e θ,0 |x × e r,0 | 2 , (A 9a) B 0 (x; x 0 ) = f θ 4πη R -r r 0 + x • e r,0 r
x • e φ,0 |x × e r,0 | 2 .

(A 9b)

Here the notation e r,0 = e r (x 0 ) was introduced for brevity. The superscript 0 denotes a free-space potential, (x; x 0 ) denotes the value of the potential at position x due to the singularity at x 0 , r = |x| and R = |x -x 0 | is the distance between the two points. For clarity, we use two ways of expressing coordinate dependence, (x; x 0 ) and (r, θ, φ; r 0 , θ 0 ), interchangeably as needed. [START_REF] Padmavathi | Motion inside a liquid sphere: internal singularities[END_REF] proved a theorem that makes it possible to construct potentials that satisfy a certain set of boundary conditions on the unit sphere from the free-space Green functions Eq. (A 9). For the free-slip boundary conditions Eq. (A 7) and Eq. (A 8) and potentials {A 0 , B 0 } that are singular only within the ball r = 1, the solution for r < 1 is given by

A(r) = A 0 (r) -rA 0 1 r (A 10a) B(r) = B 0 (r) + 1 r B 0 1 r -3r 1 r ∞ ξB 0 (ξ) dξ, ( A 10b)
where we have omitted the dependence on the variables {θ, φ; r 0 , θ 0 } for brevity.

In order for the integral to be well-defined, it is necessary that B 0 (r) decays faster than r -2 as r → ∞. However, as can be seen from an expansion of Eq. (A 9b), this is not the case for the transverse Stokeslet. In fact, the Green function for a single transverse force does not exist, because it would exert a net torque on the fluid that cannot be balanced due to the free-slip boundary condition on the sphere surface. Following [START_REF] Padmavathi | Motion inside a liquid sphere: internal singularities[END_REF], we resolve this problem by adding a rotlet singularity with strength -f θ r 0 e φ,0 at the centre of the sphere to balance the torque associated with the transverse Stokeslet. The potentials for the rotlet are

A 0 (x; x 0 ) = 0, (A 11a) B 0 (x; x 0 ) = - f θ r 0 8πη x • e φ,0 r 3 . (A 11b)
The combined potentials for a transverse Green function are then

A 0 (x; x 0 ) = f θ 8πη x • e r,0 -r 0 r 0 (R -r) + r -r 0 x • e r,0 r x • e θ,0 |x × e r,0 | 2 , (A 12a) B 0 (x; x 0 ) = f θ 4πη R -r r 0 + x • e r,0 r -r 0 |x × e r,0 | 2 2r 3 x • e φ,0 |x × e r,0 | 2 (A 12b)
The potentials in Eq. (A 12) can now be injected into Eq. (A 10) to obtain a solution that satisfies the boundary conditions at r = 1, and that behaves asymptotically like a Stokeslet near x 0 and like a rotlet near 0. In order to employ this Green function we additionally need to demonstrate that the non-local rotlet singularity does not cause a problem in the boundary-element formalism, which we do in Appendix B.

With expressions for the potentials A and B in hand, we can now use Eq. (A 5) to determine the velocities, pressure, and stress tensor associated with a unit transverse force. In order to keep the algebra manageable, we consider a 'standard' unit point force located on the x-axis at (r 0 , 0, 0) and pointing in the y-direction, and use rotations to generalise this to arbitrary positions and orientations. Let (u s x , u s y , u s z ) be the Cartesian components of the velocity due to this point force, and let p s be the pressure. Furthermore, define the auxiliary quantities c 1 = r 2 + r 2 0 -2r 0 x, c 2 = 1 + r 2 r 2 0 -2r 0 x. (A 13) a Stokeslet as usual, while as x → 0 we have

G ij ∼ Π jm ε iml x l r 3 , T ijk ∼ -3Π jm ε iml x k x l + ε kml x i x l r 5 , ( B 1)
where Π jm (x 0 ) is a projection operator that ensures that the rotlet is only present for transverse forcing. These diverge as r -2 and r -3 respectively. Hence in order to satisfy the divergence theorem, Eq. (2.3.7) in [START_REF] Pozrikidis | Boundary Integral and Singularity Methods for Linearized Viscous Flow[END_REF] needs to be augmented by an additional integral over the surface S 0 of a small ball surrounding the origin. As a result we have → 4πηω m (0), (B 4) since S0 n l dΩ = 0. Here ω = ∇×u is the vorticity of the flow. However, since the setup and flows we consider in this paper are axisymmetric, the vorticity at the origin is zero by construction. In conclusion, Eq. (B 2) reduces to the standard boundary integral equation and therefore we do not need to worry about any numerical consequences due to the additional singularity in the Green function. Physically this is because there is no net torque on the sphere for axisymmetric motion, and hence the singularities at the origin due to locally non-zero transverse forces cancel out.

Appendix C. Test problem: Concentric viscous spheres

Our goal is to develop a simple analytical flow solution that can be used to verify the correctness of our boundary-element code. We consider the flow due to a buoyant fluid sphere of radius R 1 rising/sinking inside a concentric sphere of radius R 0 with a free-slip surface containing fluid with a different viscosity. The fluid between the two spherical surfaces has viscosity η 0 and density ρ, and the fluid inside the inner sphere has viscosity η 1 ≡ γη 0 and density ρ + δρ. Gravity is directed along the Cartesian x-axis, i.e. it is not radial as it is in the Earth. While [START_REF] Happel | Low Reynolds Number Hydrodynamics[END_REF] solve a number of similar problems, the case of a free-slip outer sphere is not among them.

In the equations discussed below, lengths are non-dimensionalised by R 0 and velocities by gδρR 2 0 /η 0 , where g is the gravitational acceleration. The flow in both fluids can be described by poloidal scalars Φ i (r, θ), where i = 0 or 1 is the fluid index. The velocity in each fluid is defined in terms of the poloidal scalars as

u i = - e θ r ∂ 2 rθ (rΦ i ) + e r r B 2 Φ i , B 2 = 1 sin θ ∂ θ (sin θ ∂ θ ) (C 1)
where e θ and e r are unit vectors in the indicated directions. The poloidal scalars satisfy the biharmonic equation

∇ 4 Φ i = 0. (C 2)
Previous experience [START_REF] Happel | Low Reynolds Number Hydrodynamics[END_REF] shows that Φ i ∝ cos θ for this problem. The corresponding general solution is

Φ i = cos θ a i r + b i r -2 + c i r 3 + d i (C 3)
The pressure associated with Eq. (C 3) is

p i = -2η i (10c i r + d i r -2 ) cos θ (C 4)
where η 0 = 1 and η 1 = γ. To determine the eight constants a 0 -d 1 , we first require the solution to be finite at r = 0, which implies b 1 = d 1 = 0. Next we apply the free-slip conditions Φ 1 (1, θ) = ∂ 2 rr Φ 1 (1, θ) = 0 on the outer surface r = 1. Finally, we apply four matching conditions on the normal and tangential components of the velocity and stress at the interface r = β. The inhomogeneous condition that drives the flow is the matching condition on the normal stress σ rr , which has the form where U is dimensional. To verify the correctness of Eq. (C 10), we consider the limit β → 0 corresponding to a very small inner sphere, and re-dimensionalise U

Figure 1 :

 1 Figure 1: Definition sketch of the model. Right portion (not to scale): an axisymmetric viscous shell with thickness h, density ρ + δρ and viscosity η1 is immersed in a spherical viscous planet with radius R0, density ρ and viscosity η0. The outer surface of the planet is impermeable and free of shear traction. The shell comprises a spherical 'plate' θ θt and a downward-dipping 'slab' θt θ θ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 ≡ R0 -h/2 -d. The gravitational vector 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 R0 = 6370 km.

  φ + Ĝzr sin φ Ĝxθ Ĝyθ cos φ + Ĝzθ sin φ dφ. (4.8) Note that only the radial Green function is necessary to calculate Mrx and Mrσ , and hence the radial component Ŝ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 Dα (x 0 ) = C u β (x) Qαβγ (x, x 0 ) -u β (x 0 ) Pαβγ (x, x 0 ) n γ (x) dl(x), (4.9) where the indices β and γ take cylindrical values {x, σ} and the tensors P and Q φ + Tzαx sin φ Txαy Tyαy cos φ + Tzαy sin φ dφTyαy cos 2 φ + Tzαz sin 2 φ + Tyαz sin 2φ dφ. (4.10c) Again, only the radial (transverse) Green function is needed to calculate the radial (transverse) component of the double-layer potential at x 0 .

Figure 2 :

 2 Figure 2: Dynamics of an illustrative instantaneous subduction solution with θt = 30 • , θs = 36 • , h/R0 = 0.0157, d/h = 0.3, ϕs = 45 • and γ = 100. (a) Geometry of the shell with midsurface velocity vectors. The longest arrow has magnitude 0.299h 2 gδρ/η0. (b) Radial velocity ur. The label "fb" indicates a region of upward radial velocity corresponding to flexural bulging seaward of the trench. (c) Bending rate K and (d) stretching rate Ė 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.

  the normalised sinking speed V (St)/V Stokes predicted by the two-dimensional Cartesian BEM code of Ribe-

Figure 3 :

 3 Figure 3: Bending rate K as a function of arclength for a shell with θt = 30 • , θs = 36 • , h/R0 = 0.0157, d/h = 0.3, ϕs = 45 • and γ = 10 2 (red), 10 3 (green) and 10 4 (blue). The definition of the bending length l b is shown for γ = 10 2 .

  as functions of St and Σ for three values of the arcwise diameter D of the plate: (a) 2224 km (θ t = 10 • ), (b) 3335 km (θ t = 15 • ) and (c) 4447 km (θ t = 20 • ). Each panel is based on 50 BEM solutions with five different values of θ s -θ t and ten different values of γ.

Figure 4 :

 4 Figure 4: Normalised sinking speed V /V Stokes as a function of the flexural stiffness St for a slab with dip ϕs = 45 • . Open circles: predictions of the axisymmetric BEM code in the flat-Earth limit θt = 90 • , with d/h = 0.3 and various combinations of values of γ 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 (2010).

Figure 5 :

 5 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 Σ, for three values of the arcwise diameter D the plate: (a) 2224 km, (b) 3335 km and (c) 4447 km. Each panel comprises 50 numerical solutions with five different values of θs -θt (labelled in part (a)) and ten viscosity ratios γ ∈ [10 2 , 10 4.25 ] increasing from left to right.

  Figure 7 shows T 2 , normalised to its flat-Earth value T 2f lat , as a function of the flexural stiffness St and the sphericity number Σ. Depending on the values of St

Figure 6 :

 6 Figure 6: Hoop stress resultant T2 as a function of arclength for the same model parameters as in figure 2.

Figure 7 :Figure 8 :

 78 Figure 7: Same as figure 5, but for the hoop stress resultant T2 normalised to its flat-Earth value T 2f lat .

Figure 9 :

 9 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: CO = Cocos, JF = Juan de Fuca, NZ = Nazca, PA = Pacific, PS = Philippine Sea. The subduction zones that we consider are indicated by green lines, and their names are shown in green.

Figure 10 :

 10 Figure 10: Effect of sphericity on the sinking speed for selected terrestrial subduction zones as a function of the assumed plate/mantle viscosity ratio γ. (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 γ for the Earth.

Figure 11 :

 11 Figure 11: Definition sketch for the image system of a radial point force f located at a field point x0 with cylindrical polar coordinates (x0, σ0, 0). The colatitude of the field point is θ0 and its spherical radius is r0. The image x IM 0

S0ε

  8πηu j (x 0 ) + D [G ij (x; x 0 )σ ik (x) -ηu i (x)T ijk (x; x 0 )] n k (x) dS(x) = -Π jm (x 0 ) S0 ε iml n l r 2 σ ik (x) + 3ηu i (x) ε iml n k n l + ε kml n i n l r 3 n k r 2 dΩ, (B 2)where n k = x k /r and we recover the boundary integral formalism if the RHS vanishes. The first term in the integral on the RHS tends smoothly toσ ik (0) S0 ε iml n l n k dΩ ∝ σ ik ε imk = 0, (B 3)by symmetry of the stress tensor. Meanwhile, the second term in the integrand diverges, so we need to Taylor-expand u i (x) about the origin and use the identity ε kml n k n l = iml n l n b ∂ b u i (0) dΩ + O (r)

  (σ rr ) 0 (β, θ) -(σ rr ) 1 (β, θ) = x(β, θ) ≡ β cos θ. (C 5) The explicit expressions for the constants are b1 = d 1 = 0, d 0 = -a 0 = ββ 2 [-3 + β 5 (2β -3)(γ -1) -2γ + β(2 + 3γ)], the main parameter of interest is the ascent/descent speed U of the inner sphere. This is given by the radial (=vertical) velocity at the North pole (r, θ) = (β, 0) of the inner sphere, orU η 0 gδρR 2 0 = β 2 [2β 6 (γ -1) -β 5 (3γ -2) -2(γ + 1) + β(3γ + 2)] 3 [2 + 2β 5 (γ -1) + 3γ] ≡ G(β, γ) (C 10)

Declaration of Interests. The authors report no conflict of interest.Author ORCID. A. Chamolly, https://orcid.org/0000-0002-2383-9314; N. M. Ribe, https://orcid.org/0000-0003-0584-888X

Appendix A. Green functions

In this Appendix we derive the Green functions for the velocity and stress due to a point force inside a fluid sphere of unit radius with a constant viscosity η and a free-slip outer boundary condition on its surface. We consider separately the cases of a radially-directed force and a transverse force with no radial component. The application of the Green functions in our boundary element method is described in §4.

A.1. Radial point force

The solution for this case is surprisingly simple, and can be obtained using the method of images. Figure 11 shows a sketch of the geometry. Let x 0 be the observation or field point with colatitude θ and distance r 0 < 1 from the centre of the sphere. The image of x 0 outside the sphere is x IM 0 ≡ r -1 0 e r (x 0 ). Finally, x is another point with the same colatitude θ and radius r 0 as x 0 ; it will serve as the integration point when we perform an azimuthal integration of the Green function. The cylindrical polar coordinates of x 0 and x are (x 0 , σ 0 , φ 0 ) and (x, σ, φ), respectively. Without loss of generality, we set φ 0 = 0.

To satisfy the free-slip boundary conditions on the sphere's surface r = 1, it suffices to add to the direct Stokeslet at x 0 a second image Stokeslet of strength -r -1 0 f at the image point x IM 0 . This is shown for the case r 0 > 1 in Eq. (10.36) of [START_REF] Kim | Microhydrodynamics: Principles and Selected Applications[END_REF], but continues to hold true for an internal singularity due to the symmetry of the solution. The radial Green function Ĝir is therefore Ĝir (x;

where e r (x 0 ) = {cos θ 0 , sin θ 0 , 0} and

is the Oseen tensor. By analogy, the Green function for the stress is

where

Given the expressions Eqs. (A 1) and (A 3) for Ĝir and Tirk we must now evaluate the integral expressions for the matrix elements Mαβ , Pαβγ and Qαβγ . To do this, we write the components of Ĝ and T in terms of cylindrical coordinates by setting x -

In the foregoing expressions, the angle φ 0 has been retained for clarity, even though φ 0 = 0 because we are working in the x-y plane by hypothesis. While it is possible to write down closed-form expressions for Mαβ , Pαβγ and Qαβγ in terms of complete elliptic integrals, the results are rather complicated. Accordingly, we chose instead to evaluate the integrals numerically. This also provides for greater consistency with the case of a transverse point force (next subsection), for which one has no choice but to perform the azimuthal integrals numerically.

Then we find

)

The corresponding expressions for the components of the strain-rate tensor e s ij are not written down here due to their length, but are readily determined by differentiating Eqs. (A 14)-(A 16). The stress tensor is then σ s ij = -p s δ ij + 2ηe s ij as usual.

Finally, for a general transverse force located at position (x, y, z) = (r 0 cos θ 0 , r 0 sin θ 0 , 0), we obtain the desired tensor components Ĝiθ and Tiθk as

where

is a rotation matrix. In Eqs. (A 18) and (A 19) the quantities u s i and σ s ij are evaluated at R • x.

Appendix B. Justification of the boundary element equation with a non-local singularity

In this section we derive the boundary integral equation for our numerical scheme. Since this is largely a standard calculation we focus here on the effect of the rotlet singularity at the origin and refer the reader to [START_REF] Pozrikidis | Boundary Integral and Singularity Methods for Linearized Viscous Flow[END_REF] for details.

We define the velocity and stress Green functions G ij (x; x 0 ) and T ijk (x; x 0 ) in analogy with Pozrikidis's notation. As x → x 0 these behave asymptotically like using the radius R 1 of that sphere. We thereby find

The expression Eq. (C 11) agrees exactly with the Stokes-Hadamard-Rybczynski solution for a fluid sphere moving in an infinite fluid with a different viscosity.

Appendix D. Calculations for Pacific subduction zones

The first step in the calculations is to determine the angular radius

of a spherical cap having the same area A as the plate in question, using the plate areas from table 1 of [START_REF] Bird | An updated digital model of plate boundaries[END_REF]. The resulting values of θ t are given in column 2 of table 2.

To illustrate the subsequent steps, we begin with the Cocos plate subducting beneath Central America. We used the five trench-normal transects MEX6 and COST1-COST4 from table 1 of [START_REF] Lallemand | On the relationships between slab dip, backarc stress, upper plate absolute motion, and crustal nature in subduction zones[END_REF], which bracket a segment of the trench about 880 km long. The average values of the slab length and the deep slab dip for these transects are l = 550 km and ϕ s = 59 • , respectively. The average age of the lithosphere at the trench is τ = 21.4 Ma. The lithospheric thickness h corresponding to this value of τ can be obtained from the standard half-space cooling temperature profile

where T 0 = 0 • C, ∆T = 1325 • C, κ = 8×10 -7 m 2 s -1 is the thermal diffusivity, and z is the depth from the upper surface of the lithosphere. Defining h as the depth to the 1200 • C isotherm, we find h = 2 √ κτ erf -1 (0.9057) = 2.36 √ κτ = 55 km. The dimensionless slab length is therefore l/h = 10.0, whence our slab geometry model (2.2) implies θ s = 12.2 • given θ t = 8.65 • . The values of τ , h, l, ϕ s and θ s are given in columns 4-8 of table 2.

The calculations for the remaining subduction zones proceed similarly. In each case, we used 3 to 5 neighbouring transects from table 1 of [START_REF] Lallemand | On the relationships between slab dip, backarc stress, upper plate absolute motion, and crustal nature in subduction zones[END_REF]. These transects are TONG1-TONG4 for Tonga; NMAR1-NMAR4 for Marianas; NCHI3-NCHI6 for Chile; RYU1-RYU5 for Ryukyu; and CASC2-CASC4 for Cascadia. For Tonga and Marianas, h cannot be calculated using the halfspace cooling model Eq. (D 2), which ceases to be valid at τ ≈ 70 Ma [START_REF] Parsons | An analysis of the variation of ocean floor bathymetry and heat flow with age[END_REF]. For greater ages the thickness of the lithosphere tends to a constant asymptotic value h ≈ 100 km, and we therefore used this value for Tonga and Marianas.