We present a novel approach for modeling subduction using a Multipole-accelerated Boundary Element Method (BEM). The present approach allows large-scale modeling with a reduced number of elements and scales linearly with the problem size. For the first time the BEM has been applied to a subduction model in a spherical planet with an upper-lower mantle discontinuity, in conjunction with a free-surface mesh algorithm.