Subduction dynamics is strongly dependent on the geometry and rheology of the subducting slab and adjacent plates, as well as on the induced mantle flow driven by the evolution of tectonic configurations along subduction zones. However, these processes, and the associated plate tectonic driving forces, are difficult to study using time‐dependent 3‐dimensional computer simulations due to limitations in computing resources. We investigate these phenomena with a novel numerical approach, using BEM‐Earth, a Stokes flow solver based on the Boundary Element Method (BEM) with a Fast‐Multipole (FM) implementation. The initial BEM‐Earth model configurations self‐consistently determine the evolution of the entire lithosphere‐mantle system without imposing additional constraints in a whole‐Earth spherical setting. We find that models without an overriding plate overestimate trench retreat by 65% in a 20 m.y. model run. Also, higher viscosity overriding plates are associated with higher velocity subducting slabs, analogue to faster oceanic plates subducting beneath more rigid continental lithosphere. In our models poloidal flows dominate the coupling between the down‐going and overriding plates, with trench‐orthogonal length variations in overriding plates inducing flows at least ∼2× stronger than trench‐parallel width variations. However, deformation in the overriding plate is related to its length and width, with narrower and longer plates extending more than wider and shorter plates. Key Points Quantify the role an overriding plate has on subduction dynamics Determine subduction zone geometry and rheology influence on plate velocities Determine ratios of driving forces distributed in simplified subduction systems