In subduction zones, the subduction of oceanic plates would drive complex mantle flow with a pattern affected by several factors, including the thermal structure of the subducting slab, the slab shape, the overriding plate thickness, and the viscosity of the surrounding mantle (e.g., Funiciello et al., 2003;Jadamec & Billen, 2010;Stegman et al., 2006). For a single slab subduction system, the dominant flow regimes include 2-D corner flow in the mantle wedge and subslab domain induced by viscous coupling between the subducting slab and the surrounding mantle (e.g., Kincaid & Sacks, 1997;Van Keken et al., 2002), and 3-D toroidal flow invoked by slab rollback (e.g., Faccenda & Capitanio, 2012;Long & Silver, 2008). The flow pattern may become more complicated when considering interactions of multiple subducting slabs as revealed by recent numerical modeling studies (e.g., Di Leo et al., 2014;Holt et al., 2017;Király et al., 2018). The eastern SE Asia region represents one of the ideal places for studying dynamics and evolution of subduction zones because of its long-lived subduction history and the assembly of slab subduction at various stages (Figure 1) (e.g., Hall, 2002;Yin, 2010). In its northern part, the Philippine region is characterized by the presence of active subduction zones on its both sides. To the west, the South China Sea (SCS), Sulu Sea, and Celebes Basins are subducting beneath the Philippine archipelago along