A Possible Roll-Over Slab Geometry Under the Caroline Plate Imaged by Monte Carlo Finite-Frequency Traveltime Inversion of Teleseismic SS Phases

The shape of a subducting slab varies as a function of trench motion. Two end-members of subduction modes are geodynamically possible: roll-back mode underneath neighboring plates and roll-over mode underneath the plate itself. Whereas most of major slabs seem to roll back while the Pacific plate shows a slab piling behavior down to ∼1,000 km depth under the Mariana trench, no clear evidence of slab roll-over in nature has been reported so far. Here we show a possible roll-over slab beneath the Caroline microplate, revealed from its three-dimensional seismic velocity structure derived by analyzing teleseismic reverberating SS phases. We suggest that slab roll-over is driven by at least two factors: 1) the overall buoyancy and fragility of the Caroline microplate at the surface, induced by a thin hot mantle plume that rises from depths ≥800 km; and 2) the pushing force of the Pacific plate acting on the trailing edge of the Caroline plate.


INTRODUCTION
Seismic tomographic studies have imaged various subduction patterns (e.g., Okino et al., 1989;Fukao and Obayashi, 2013;French and Romanowicz, 2015;Hosseini et al., 2020). Whereas a slab with a retreating trench is dominant for the Pacific plate going under the North American plate in the north-eastern Japan arc (e.g., Lallemand et al., 2008;Fukao and Obayashi, 2013), a steep slab with a slab piling geometry of recumbent folds is also possible, as below the advancing Mariana trench (e.g., Widiyantoro et al., 1999;Barklage et al., 2015). These observed subduction modes have motivated the community of geodynamics to foster the interpretation of slab physics, with the aid of global and regional numerical simulation and analogue laboratory experiments. Geodynamical modeling studies, both numerical and experimental, show that the geometry of subducting slabs is closely related to the motion of converging plates (Christensen, 1996;Schellart, 2011): the roll-back mode occurs when trench moves oceanward ("retreats"), whereas the roll-over mode occurs when trench moves landward ("advances"; Replumaz et al., 2004;Manea and Gurnis, 2007) (Schellart, 2008). It is also shown that slabs may not only roll back beneath neighboring plates, as has been observed beneath the Mariana trench, but that it could also roll over below themselves (Schellart, 2008).
One of the possible mechanisms for a slab to roll over is by a large pushing force along the trailing edge of a subducting oceanic plate (Schellart, 2005;Stegman et al., 2010a). This mechanism may not apply to the major subducting plates on the Earth, given that the pushing force at the trailing-edge of large, old oceanic plates is normally much smaller than their gravity, which is a principal engine of plate motions (e.g., Forsyth and Uyeda, 1975). In addition, heavy, stiff, old subducting slabs may not achieve the slab roll-over mode, which is likely to occur in buoyant and highly deformable slabs (Stegman et al., 2010b). Hot plumes can contribute to such pushing forces by providing heat under oceanic plates, thereby making them more buoyant and deformable. Microplates are thus worth exploring to better understand the geodynamics of subducting slabs which may exhibit roll-over mode.
In this paper, we explore the Caroline microplate as a good candidate for slab roll-over mode. This oceanic microplate is surrounded by three major old plates and fed by a plume. We investigate the mode of subduction of this small plate by first developing and applying a new seismological imaging technique particularly adapted to the region of interest. We then propose the geodynamical model to reconcile the seismological and tectonic models.

BACKGROUND Tectonic Setting of the Caroline Plate
The Caroline plate has a complex and discussed geologic history, that includes tectonic rotations, protracted and recent subduction, former rift jumping and abortion, young spreading centers, plate boundaries of unknown geometry and kinematics, and hotspot and volcanic activity (Gaina and Muller, 2007). This microplate is roughly hexagonal (approximately 1,800 km × 1,000 km) and is trapped among three major converging plates (Figure 1): the Australian (south), Philippine (west and north-west), and Pacific (east and northeast) plates (Weissel and Anderson, 1978). The Philippine-Caroline plate boundary, located in the west, consists of convergent boundaries to the north, marked as the Yap and Palau trenches, and the Ayu trough extensional plate boundary to the south (Weissel and Anderson, 1978). The south and east boundaries of the Caroline plate evolve in an obliquely and rapidly converging plate boundary zone (Baldwin et al., 2012). Tectonic accommodation seems to occur by lithospheric shortening along the New Guinea, Manus, and Mussau trenches. The Caroline plate is underthrusting beneath the Indo-Australian and Pacific plates. To the north, its contact with the Pacific plate is unclear but may be dominated by a spreading center controlling the evolution of the Sorol trough (Weissel and Anderson, 1978). The transition between this plate FIGURE 1 | Tectonic setting of the Caroline plate. The lower-right inset shows the larger tectonic framework of the Caroline plate. The main panel depicts the nature of the main tectonic boundaries around the Caroline plate (Cooper and Taylor, 1987;Gaina and Muller, 2007;Baldwin et al., 2012). Ph, Pa, and InAu are the acronyms for Philippine, Pacific, and Indo-Australian plates, respectively, the three major plates surrounding the Caroline Plate. WCT and KT are the West Caroline and Kiilsgaard troughs, remnants of a former spreading center. CHT stands for Caroline Hotspot Track. The bathymetric shaded-relief visualization of the GEBCO 2014 grid is shown as produced by NOAA/NCEI. Bathymetric depths and major morphologies are marked with equal-depth contours in kilometers.
boundary and the Mussau trench is of complexity, potentially due to its interaction with the Caroline hotspot track (CHT in Figure 1; Gaina and Muller, 2007;Wu et al., 2016). The dating of the CHT could be the trace of the Caroline plume from the core-mantle boundary (CMB), which might have been a source for the Caroline ridge, and probably the Sorol trough.
Subducting features of the Caroline plate are observable at its southern plate boundaries, along the New Guinea and Manus trenches. Although the lack of geodetic data hinders further details, east-west trending magnetic anomalies within the Caroline plate seafloor indicate southward plate motion (Maus, 2009). A north-south volcanic ridge, called the Eauripik rise, rises from the consistently flat seafloor bathymetry at ∼4.5 km depth, while the seafloor shallows by more than 2 km around the center of the Caroline plate ( Figure 1; Weissel and Anderson, 1978). The Eauripik rise has a relatively constant width of ∼250 km in its northern and central sectors that triples in its southern sector, as the ridge becomes shallower and more distributed (Figure 1). The Eauripik rise disrupted the north-south spreading center that originally formed this oceanic microplate and defined the east and west Caroline basins during the Oligocene (Weissel and Anderson, 1978;Gaina and Muller, 2007). This rise was once thought to be a dead divergent boundary (Winterer et al., 1971). However, the magnetic exploration and borehole analyses suggest that the Eauripik rise may be related to magma that leaked through faults and cracks when the plate was created (Weissel and Anderson, 1978;Hegarty and Weissel, 1988). Another interpretation suggests a link between the Eauripik rise and the tracks of the Manus plume (Macpherson and Hall, 2001). Subsurface and seismic data leading to high-resolution imaging can resolve the mantle structure underlying this complex tectonic framework at present, and help understanding potential links between lithospheric and surface processes.

Previous Seismological Observations
Looking from the bottom of the mantle, global whole-mantle waveform inversion studies indicate the existence of a plume with a maximum size of 10°× 10° (French and Romanowicz, 2015;Garnero et al., 2016). As imaged by these studies, the Caroline plume is a vertical pile of low shear velocity (V S ) rooted at the western tip of the Pacific Large Low Shear Velocity Province. Body-wave waveform inversion studies imaged the detailed elastic and anelastic structure of the plume, suggesting that the Caroline plume is iron rich at the base of mantle (at least until ∼1,000 km above the CMB), probably due to chemical interaction with the CMB (Konishi et al., 2017;Deschamps et al., 2019;Konishi et al., 2020). It is thus interesting, from the geodynamical point of view as well, to trace this Caroline plume up to the surface and to explore its link to the Caroline ridge.
The fine structure of the upper mantle and mantle transition zone beneath the Caroline plate has not been imaged, for there is not any dense array of seismometers nor seismicity around the region except near its south-western boundary. Long-wavelength global tomography models have overcome this difficulty and, using all the seismic phases that sample any region of the globe, have successfully shown a large high-velocity anomaly in the transition zone beneath the Caroline plate (e.g., O'Neill et al., 2005;French and Romanowicz, 2015). However, the detailed structure within the region is still unrevealed. Several factors are responsible for the lack of resolution: 1) previous studies use lowfrequency contents (shortest periods of 32-128 s, i.e. larger than ∼150 km) in their inversions and thus have a poor sensitivity in resolving ∼100 km-scale structure; and 2) due to a heavy regularization, smoothing scheme and data weighting factors, to stabilize inversions for a number of parameters within the whole Earth domain, it can bias the obtained results and broaden the anomaly images (e.g., Bozdag et al., 2016;Marjanović et al., 2017). Global tomography can also prevent from imaging sharp interfaces that can be observable from forward modeling (e.g., Ni et al., 2005). This results in the difficulty to assess the error for each single point in global tomographic models. In order to overcome these limitations, we develop a Monte Carlo (MC) inversion using SS-S differential traveltimes and waveforms, specially designed for the purpose of exploration of the interior of the upper mantle and the transition zone beneath the Caroline plate. We apply it to a large dataset of waveforms that sample the region of interest.
Given the distribution of teleseismic earthquakes (Supplementary Figure S1), we provide the first threedimensional (3D) mantle structure to reveal slab roll-over geometry beneath the Caroline plate and understand subduction dynamics. On the basis of our model, we aim to understand the structure and evolution of the Caroline plate, which are affected by several major plates that collide and subduct one beneath another. We will further focus on the following questions with the aid of seismology and geodynamics: 1) which mechanism is supplying the materials at the Caroline ridge and what is the link between the hot iron-rich Caroline plume at the mantle (e.g., O'Neill et al., 2005;Konishi et al., 2017;Deschamps et al., 2019;Konishi et al., 2020)? and 2) how has the Eauripik rise been formed at the center of the Caroline plate?

Monte Carlo Inversions Using SS-S-phase Double-Difference Traveltime and Waveform Similarity
Teleseismic waveform data that sample only inside the upper mantle and the mantle transition zone beneath the Caroline plate will not provide enough information on its structure, given the small number of seismic stations and earthquakes inside the region. We thus specifically chose an SS seismic phase for retrieving V S structure, for it reverberates at the surface of the plate once between the source and the receiver. As the direct S and SS phases sample similar regions at both the source and receiver sides, the differential traveltime between their arrivals is mainly due to the bouncing point of the SS phase. We aim to obtain preferred 1D model(s) for the bouncing point region by using a MC approach which retrieves as much robust information as possible from the data and avoids biasing the inversion with a priori information, such as initial models.
We collected a dataset consisting of 2,540 transversecomponent waveforms from 482 teleseismic earthquakes whose SS phase bounced within and around the Caroline plate (Supplementary Figure S1). The collected events have Mw 5.5-7.0 and focal depths of <75 km, and occurred in 2001-2015. We filtered the observed waveforms at 20-100 s, and the data were quality checked. We retained only waveforms with a signal-to-noise ratio exceeding 2.5 for S and 1.5 for SS, and divided the dataset based on bouncing points of the SS phase. We computed traveltime kernels for S and SS arrival times of the data and their differences using DSM Kernel Suite (Fuji et al., 2012;Fuji et al., 2016). For epicentral distances of 50-90°, the sensitivity of SS phase is focused inside a cylindrical column of a radius of 5°beneath the bouncing point. The kernels show that SS phases exhibit significant sensitivity to structure, indicating a possibility of extracting information on localized inner structure of SS bouncing point region (Supplementary Table S2). We note that we cannot fully eliminate influences to SS traveltimes from the source-side, receiver-side and lower mantle heterogeneities even with the data covering a wide range of azimuths (Supplementary Figure S1). We thus test statistical robustness of our inversion a posteriori as a function of the number of models chosen, that are represented as covariance of our inverted models. We further note that we did not use any partial derivatives or kernels during the inversions, although both can be the indices for the determination of an adequate binning of sub-regions and for that of frequency ranges.
We then generated synthetic seismograms for a family of 1,000 1D models of V S with a ±4% anomaly in the depth range of 0-850 km for the Preliminary Reference Earth Model (PREM; Dziewonski and Anderson, 1981), allowing discontinuities anywhere. We fixed the number of layers to be 9, but the velocity (and its gradient) of each layer can be changed randomly (Supplementary Figure S2). We computed the synthetic seismograms using the 1D direct solution method (Geller and Ohminato, 1994;Geller and Takeuchi, 1995). The MC inversion was aimed at finding a 1D model(s) for each subregion by minimizing the misfit function: where T obs/syn phase denotes the traveltimes of each phase for the observed and synthetic phase waveforms (for a model m). The first L1-norm term corresponds to the difference of SS-S differential traveltimes between the synthetic and observed data, denoted here as double-difference traveltime (DDTT). We measured the differential traveltimes for S and SS phases by cross-correlating the observed and synthetic seismograms. CC denotes the cross-correlation coefficients, α is a weighting factor that we set empirically to 5.0, and N d is the number of waveforms.
Looking at various epicentral distance seismograms has the advantage that their sensitivity of SS phases with respect to each depth differs from one to another, and thus results in a quasi-orthogonal sensitivity, when using a large number of waveforms (Konishi et al., 2014). Furthermore, our MC method allows us to not depend on a priori information, or an initial model, resulting in a global minimum search of the misfit function (Eq. 1). As shown in the Results section, we also evaluate the covariances for our preferred models, thereby providing an index of confidence level of each model parameter.

Geodynamic Simulation
We conducted numerical simulations to constrain controlling factors in the geometry of a subducting microplate. We considered a 2D trench-normal cross section, as trenchnormal motion between converging plates is the primary factor in the geometry of subducting slabs (e.g., Christensen, 1996;Torii and Yoshioka, 2007;Stegman et al., 2010a;Stegman et al., 2010b;Schellart, 2011;Cizkova and Bina, 2015). We set up a 2D model of a subducting slab with an infinitely large width, given that 3D simulations of subduction conclude that slab behavior is sustained when their width is >1200 km, which is the case for the Caroline plate (Stegman et al., 2010a). Our model includes an overlying plate, a subducting microplate, and a pushing plate in a whole-mantle domain (Supplementary Figure S3). Each plate is decoupled by fracture zones characterized by weak strength (e.g., Manea et al., 2014). In addition, low-viscosity weak zones were set at the top-left and top-right corners of our model domain to allow horizontal plate motions. Each plate consists of overlying crustal and underlying mantle rocks with their proper densities and heat generation rates (e.g., Turcotte and Schubert, 2002). The model includes two discontinuities associated with the olivine/wadsleyite and ringwoodite/bridgmanite phase changes: the subducting microplate acquires negative buoyancy at the 410 km depth and positive buoyancy at the 660 km depth, resulting in slab deformation.
We used basic equations and numerical solvers of previous studies (Nakao et al., 2016;Nakao et al., 2018): mantle flows were derived from a Stokes equation for an incompressible viscous continuum with an infinite Prandtl number; and mantle temperature was solved using the extended Boussinesq approximation (see Supplementary Material text). Note that, contrary to the previous studies, we excluded water transport of the simulation for simplification. In order to solve mantle flows along four model boundaries, free slip conditions are imposed, except for the top of the pushing plate. Along the pushing plate surface, a landward constant velocity (V push ) is imposed instead. We set a constant temperature of 0°C as boundary conditions along the top boundary, whereas the insulating condition is imposed at the bottom, right and left boundaries (Supplementary Figure S3). We set initial conditions as follows: 1) the overlying plate is 10 Myr old and 3,500 km in length; 2) the subducting microplate is 10 or 50 Myr old and 2,000 km in length; 3) the pushing plate is 100 Myr old and 3,700 km in length; 4) the fracture zone is placed between the overlying and the subducting plates; and 5) another fracture zone is assigned between the subducting and the pushing plates (Supplementary Figure S3). Lastly, in these models, we varied the pushing velocity V push (3, 6, or 9 cm/yr) and the initial Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 593947 thermal age of the subducting microplate (10 or 50 Myr) as variable parameters to examine subduction modes.

RESULTS
Individual DDTT values show that SS phases arrive much faster than PREM prediction and S phases arrive as predicted by PREM (Supplementary Figure S4A). The averaged DDTT over the whole region is 7.6 s, implying the presence of a fast V S anomaly beneath the Caroline plate, which is consistent with global studies (e.g., O'Neill et al., 2005;French and Romanowicz, 2015). DDTT variations as a function of SS bouncing points are relatively small (<2 s) at the southern and western plate boundaries (Supplementary Figure S4B). Our MC inversion outputs preferred depth-dependent V S models for each sub-region, defined by the reverberating points of SS phases. Figure 2 shows a final seismic velocity model compiled from all preferred 1D models. We obtain error estimates by means of the MC inversion scheme that directly compares observed and synthetic data generated for 1,000 models. The final model shows high-velocity anomalies throughout the upper mantle and the transition zone, which are, again, consistent with the global models. However, more detailed small-scale features are visible in our reconstructed 3D model due to the shorter-period contents and our optimized inversion scheme that maximized resolving power and particularly targeted the inner structure beneath the SS bouncing points.
Faster velocity features (1.6-2.2% faster than the PREM) are clear, especially around the 200-750 km depths (Figure 2A, Supplementary Figures S5-S7) and mostly concentrated below the southern plate boundaries of Caroline (Figure 2A, Supplementary Figures S5-S8). Below the northern plate boundaries, fast V S features observed around the 200-850 km depths can be also recovered from the best 10 and 30 models (Supplementary Figures S9A, S10A, respectively). On the other hand, a smaller-scale low-V S structure (0.5-1.1% slower than PREM) is observed from bottom to top of our model domain (0-850 km depths). Low V S features are more prominent at shallower depths (0-200 km) and are dominant merely below the Eauripik rise at 50 km depth ( Figure 2B, Supplementary  Figures S5-S8, S9B, S10B). The standard deviations for both 10 and 30 preferred models are <1.5% beneath the Caroline plate, supporting the confidence level of our model ( Supplementary  Figures S11, S12).
Our 2D geodynamic simulation demonstrates that the geometry of a subducting oceanic plate is critically dependent on the magnitude of landward velocity (V push ) and the age of the slab (Figure 3). Subduction of a 50 Myr-old oceanic plate can lead FIGURE 2 | Three-dimensional S-wave velocity (V S ) structure beneath the Caroline plate down to 850 km depth with interpretation (arrows). The green perimeter at the surface (0 km depth) indicates the Caroline plate boundaries. (A) The region in blue has V S 2.2 and 1.6% faster than PREM (Dziewonski and Anderson, 1981). (B) The region in red has V S 1.1 and 0.5% slower than PREM. Our model suggests the presence of subducting and stagnant slabs (blue regions) from the neighboring plate boundaries, and the plume (red regions) emerging from the base of the model domain (850 km depth) and extending to the surface. See Supplementary Figure  S5  to a variety of slab geometries, depending on the V push in the range for actual oceanic plate motions (Lallemand et al., 2008). A small V push (3 cm/year) results in a slab-draping geometry with a horizontal slab segment in the transition zone ( Figure 3A), similar to the Pacific slab geometry in north-eastern Japan (e.g., Huang and Zhao, 2006). An intermediate V push (6 cm/ year) develops a draping geometry with recumbent folds ( Figure 3B), similar to that in the southern Izu-Bonin arc (e.g., Huang and Zhao, 2006). A large V push (9 cm/year) develops a slab piling geometry of recumbent folds at the transition zone ( Figure 3C), similar to the one in the Mariana arc (Widiyantoro et al., 1999). The transition between the above-mentioned slab geometries is different for the case of a young, hot oceanic plate. Subduction of a 10 Myr-old plate with a fast velocity (V push 3 cm/year) exhibits a draping geometry with recumbent folds ( Figure 3D). With intermediate to fast velocities (V push 6 and 9 cm/year), the 10 Myr-old plate exhibits a roll-over geometry with a horizontal slab segment in the transition zone, resulting in trench advance ( Figures  3E,F). In this case, the oceanic plate is continuously supplied to the converging plate boundary due to the large pushing velocity, whereas the convergence velocity slows down by the large thermal buoyancy of the subducting slab. Consequently, the overlying plate is pushed backward, resulting in the advancing trench and the roll-over geometry.
To sum up, our seismic model indicates possible roll-over slab subductions from the southern boundary and feasibly also from the northern boundary, as well as the presence of the plume-like low-velocity anomalies extending from (at least) 800 km depth beneath the Caroline plate. Our geodynamical simulations further demonstrate that the roll-over subduction needs a young, hot oceanic plate and relatively fast landward pushing velocity.

DISCUSSION AND CONCLUSIONS
Our seismic model, derived from MC inversion of teleseismic SS phases, shows a high V S feature in the upper mantle and transition zone beneath the Caroline oceanic microplate. We observe a nearvertically continuous high-velocity anomaly throughout the upper mantle ( Figure 2, Supplementary Figure S8, green arrow), and we consider it as a subducting slab from the southern boundaries of the Caroline plate. This "slab" seems to fall down vertically, deflect at 660 km discontinuity and advance backwards beneath the plate itself in the transition zone. Meanwhile, geodynamical modeling predicts that a young, hot oceanic plate with a fast pushing velocity can roll over, and such a geometry is inferred from our seismic model. On this basis, our model could provide a rare observation of a case of roll-over subduction due to the particular plate geometry surrounding the Caroline plate. The geometry of subducted slabs depends critically on the partitioning of the convergent velocity at the surface into its subducting plate motion component and trench migration component (Schellart, 2011), as confirmed here in this study. Hence, the pushing force due to the Pacific plate motion may have been enhanced, resulting in roll-over plate geometry. The roll-over mode is uncommon in actual subduction zones but is predicted to appear where 1) a subducting plate is very hot; and 2) the trailing edge of the subducting plate is vigorously pushed by another oceanic plate.
Another notable feature in our seismic model is the slowvelocity material that ascends from the bottom of the model space and detours around the slabs, reaching the surface ( Figure 2B, magenta arrow). This slow-velocity feature might be the upper end of the Caroline plume, coming either directly from the CMB or rooted at some shallower depth. This plume may be significantly less rigid than the subducting slabs and thus rise upwards through the upper mantle. Its close proximity to the Eauripik rise may indicate that the rise provides a magma   Steinberger and O'Connell, 1998). Based on our seismic model ( Figure 2) and the geodynamic simulation results (Figure 3), we developed a schematic interpretation (Figure 4, Supplementary Figure S13). The Caroline plume rises upward from the deep mantle, detouring around the cold and rigid bodies and perhaps feeding the Caroline ridge (Supplementary Figure S13, magenta arrow). The Caroline ridge, that spreads at higher rates (Muller et al., 2008), causes the plate to move southward. The microplate is probably pushed southward by Pacific plate and fed by the mantle plumes, and subducts at its southern boundaries, possibly resulting in a roll-over subduction mode beneath the microplate itself. However, the shape of the slab depends on the tectonic history and there are other possible scenarios that can create the imaged slab shape (see Supplementary Material text).
To date, the active effect of mantle plumes on the subducting slabs (e.g., Chang et al., 2016) and interactions with the mantle and neighboring slabs (e.g., Kiraly et al., 2018) remain poorly understood. Thus, the link between the plume from the CMB and the complex features of the upper mantle structure beneath the Caroline plate remain unclear. Hence, the ocean-bottom geophysical data on the plate (e.g., "Pacific Array"; Forsyth and Detrick, 2003;Kawakatsu et al., 2016) will help to better probe the upper-mantle interaction between plumes and slabs. As we showed in our models, an interplay between the upwelling plume and the oceanic microplates can have features uncommon to major plates, whose dynamics are not yet well constrained.

DATA AVAILABILITY STATEMENT
Seismic waveform data can be obtained from the IRIS Data Management Center at www.iris.edu. See Supplementary Tables S3, S4 for data we used in this study.

AUTHOR CONTRIBUTIONS
NF and S-ML organized binational workshops between the Université de Paris and Seoul National University on the 13th of November 2015, in Paris and on the 3rd of May 2016, in Seoul, to set up the project. NF developed and designed the research, and HJ carried out the data collection and analysis under the supervision of NF, YHK, and S-ML. AN conducted geodynamic modelling. HJ, NF, and YHK prepared the first draft. DF-B interpreted the tectonic setting of the region, did Figure 1, and improved the overall quality of manuscript, and AS provided the schematic interpretation of our model. S-ML provided thorough knowledge about the Caroline plate and encouraged all of the authors to investigate the detailed mantle structure beneath the plate. KK and NF maintained the seismological numerical tools for the analysis conducted in this work.   that the Caroline plate is moving southward, subducting under the Australian plate. The ridge push force is much larger than the trench migration velocity from the Australian plate, due to a relatively large volume of plume supplied from the base. This can lead to an advancing subduction regime and a roll-over slab geometry, as observed in our seismic velocity model ( Figure 2). The Caroline plume upwelling from the core-mantle boundary or deep mantle is hot and soft, and thus takes a detouring path that avoids the subducting slabs within the region. See Supplementary Figure S13 for interpretation.
Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 9 | Article 593947 Nakakuki for providing his numerical code. Numerical simulations were conducted by DA systems of the Japan Agency for Marine-Earth Science and Technology (JAMSTEC). Numerical computations were performed on the S-CAPAD platform, IPGP, France.