We present the first application in geodynamics of a (Fast Multipole) Accelerated Boundary Element Method (Accelerated-BEM) for Stokes flow. The approach offers the advantages of a reduced number of computational elements and linear scaling with the problem size. We show that this numerical method can be fruitfully applied for the simulation of several geodynamic systems at the planetary scale in spherical coordinates, and we suggest a general approach for modeling combined mantle convection and plate tectonics. The first part of the paper is devoted to the technical exposition of the new approach, while the second part focuses on the effect played by Earth curvature on the subduction of a very wide oceanic lithosphere (W = 6,000 km and W = 9,000 km), comparing the effects of two different planetary radii (ER = 6,371 km, 2ER = 2 9 6,371 km), corresponding to an ‘‘Earth-like’’ model (ER) and to a ‘‘flat Earth’’ one
(2ER). The results show a distinct difference between the two models: while the slab on a ‘‘flat Earth’’ shows a slight
undulation, the same subducting plate on the ‘‘Earth-like’’ setting presents a dual behavior characterized by concave
curvature at the edges and by a folding with wavelength of the order of magnitude of 1,000 km at the center of the
slab.