TY - JOUR
AU - Kaluza,, Owen
AB - SUMMARY In this paper, we present 3-D numerical simulations in which a compositional mantle plume rises through a shallow mantle layer at the rear of a retreating slab. The slab–plume buoyancy flux ratio Bs/Bp is varied from 3 to 30 over nine slab–plume simulations. The plume causes an overall decrease of the slab retreat rate from 15 to 7 % in the corresponding range |$B_s/B_p=[3,30]$|⁠. The retreat rate decrease occurs in two stages: the first decrease occurs remotely when plume and slab are hundreds of kilometres apart; the second decrease is linked to the slab–plume-head impact. Continuous tracking of key positions along the plume head and conduit together with velocity profiling further shows a very close interplay between the conduit and plume head dynamics. In particular, the combination of the slab– and plume(head)–induced flows at the rear of the slab increases the advection and tilt of the conduit and causes its flaring with height in the direction parallel to the trench. As a result, the conduit source slowly drifts away from the slab by hundreds of kilometres and flares by one and a half times its original radius over 30 Myr of plume head spreading. The conduit tilt and flaring result in an increasingly unbalanced azimuthal distribution of the incoming plume flux from the feeding conduit into the head. These changes in the feeding conditions coupled with the mantle flow at the level of the plume head lead to the asymmetric spreading of plume material beneath the plate in the preferential direction that is parallel to the trench. Upon its arrival against the slab, it can be a front of buoyant material set to subduct along the slab width that has widened up to two and a half times its initial dimension. In nature, it is expected to extend from hundreds to thousands of kilometres depending on the slab–plume buoyancy flux ratio. To our knowledge, this study is the first to highlight the inter-relationships between a plume head and its feeding conditions in the plume source region. Mantle processes, Numerical modelling, Dynamics: convection currents, and mantle plumes, Subduction zone processes 1 INTRODUCTION The geological record is marked by events during which mantle plumes interacted with subduction zones. The interaction has raised interests because it is often associated with the genesis of valuable deposits (e.g. Tassara et al. 2017). Notable examples include plume modified orogenies in eastern Precambrian Australia (Betts et al. 2009), in the western United States (Brendan Murphy et al. 1998), in the Abitibi–Wawa province (Wyman et al. 2002) or in the southern Central Andes (Gianni et al. 2017). In these cases, the interaction involved a geodynamic setting where a plume head is beneath a subducting plate in retreating motion, which eventually leads to its impact with the slab. Several studies numerically simulating such interaction showed that a plume head could change the geometry of a trench, modify the trench motion rate, even switching the trench retreat motion to an advancing motion, and last but not least, create a window/tear within the slab allowing mantle plume material to reach the mantle wedge and the overriding plate region (Betts et al. 2012; Lee & Lim 2014; Chang et al. 2016; Leonard & Liu 2016). However, so far, much less attention has been given to the plume feeding conduit: Had the conduit detached from the plume head at the time of interaction or was the plume head still connected to its conduit? The two scenarios in effect imply different dynamics: in the first case, the plume head constitutes a fixed-volume spreading under its own buoyancy and the ambient mantle flow (e.g. Betts et al. 2012; Lee & Lim 2014). In the second case, the plume conduit releases a flux of buoyant material into the head, which then spreads under its buoyancy and the ambient mantle flow (e.g. Chang et al. 2016). The rates of plume spreading in the two cases are expected to differ significantly. For instance, in the case of an axisymmetric viscous gravity current, the growth of the radius varies as tα with a propagation that is three to four times slower under the fixed volume condition than the fixed flux condition (Huppert 1982; Lister & Kerr 1989). Besides, as the conduit itself is also under the influence of the ambient mantle flow, the conditions of release may change with time, and in return modify plume head spreading. Such details of the supply of material by the plume conduit have been ignored in previous studies of mantle plume interaction with the lithosphere (e.g. Olson 1990; Ribe & Christensen 1994; Sleep 1996; Chang et al. 2016; Leonard & Liu 2016). Previous laboratory experiments provided the first insights into the dynamics of a head-to-toe plume overridden by a retreating slab. In particular, a decrease of the retreat rate was found for subduction in the presence of a plume of high buoyancy flux while the conduit was found to be deflected and ultimately flattened (Mériaux et al. 2015, 2016). The full extent of the inter-relationship between the plume head spreading and feeding conduit dynamics was however not addressed. The motivation of this study is to fill this gap by using numerical modelling as the method of choice because the approach simultaneously gives access to quantitative metrics and the velocity flow field in the domain, especially at the rear of the slab. 2 METHODS 2.1 Model setup The model is composed of a plate lying on a shallow mantle layer and fixed at its trailing edge, an ‘air’ layer overlying the plate and mantle wedge region, which does not include an overriding plate, and a compositional plume that is located at the base of the mantle and rear of the slab tip in the symmetry plane of the plate (Fig. 1). The plate of thickness d and width W is negatively buoyant relative to the mantle with density contrast Δρs and it has an iso-viscous Newtonian viscosity μs. The plume fluid density is less than that of the mantle by Δρp, and it has a Newtonian viscosity μp. Subduction is initiated by prescribing a pre-existing kink of the plate in the mantle that is the initial slab tip. The subduction evolution is then self-sustained by the negative buoyancy of the lengthening slab. Due to the shallow depth and fixed plate trailing edge, the subduction dynamics is mostly characterized by a retreat of plate hinge. We thus defined the slab buoyancy flux Bs by $$\begin{eqnarray*} B_s/g=\Delta \rho _s |u_s^r| W d, \end{eqnarray*}$$(1) where |$|u_s^r|$| is the slab steady retreat rate and g is the gravitational acceleration. The initiation of the plume inflow is synchronous to that of subduction. The plume inflow is characterized by a vertical Poiseuille velocity profile |$v(r)$| centered at the plume source position (xp, 0, 0) and distributed over the area |$\pi r_p^2$|⁠, where rp is the plume radius. The vertical velocity |$v(r)$| is given by $$\begin{eqnarray*} &v(r)=\frac{\Delta \rho _p g}{4\mu _p} \left( r_p^2 -r^2 \right) & \mbox{, for } r \le r_p, \end{eqnarray*}$$(2) $$\begin{eqnarray*} &v(r)=0 & \mbox{, for } r \gt r_p, \end{eqnarray*}$$(3) where |$r^2=(x-x_p)^2 + z^2$|⁠. Integration over r of eq. (3), gives the plume volume flux $$\begin{eqnarray*} Q=\frac{ \pi \Delta \rho _p g r_p^4}{8\mu _p}. \end{eqnarray*}$$(4)Q is fixed throughout the entire simulation. The plume buoyancy flux Bp is then given by $$\begin{eqnarray*} B_p/g=\Delta \rho _pQ. \end{eqnarray*}$$(5) Figure 1. Open in new tabDownload slide Generic numerical model. The two quantities yt(t) and xh(t) are two of the metrics that will be used for the quantitative estimate of the sinking of the slab tip to the bottom boundary, and the plate hinge retreat. The position xc(y, t) is the horizontal position of the centre of mass of the plume material in the plane of symmetry z = 0 and height y and the position |$z_{max}(y,t)$| is the plume material maximum extent in the z direction at height y. The positions xd(t) and xu(t) are the plume fronts downstream and upstream of the slab retreat, respectively, in the plane z = 0. The plume head maximum position beneath the plate and parallel to the trench is zl(t). Figure 1. Open in new tabDownload slide Generic numerical model. The two quantities yt(t) and xh(t) are two of the metrics that will be used for the quantitative estimate of the sinking of the slab tip to the bottom boundary, and the plate hinge retreat. The position xc(y, t) is the horizontal position of the centre of mass of the plume material in the plane of symmetry z = 0 and height y and the position |$z_{max}(y,t)$| is the plume material maximum extent in the z direction at height y. The positions xd(t) and xu(t) are the plume fronts downstream and upstream of the slab retreat, respectively, in the plane z = 0. The plume head maximum position beneath the plate and parallel to the trench is zl(t). The ratio of the slab–plume buoyancy flux Bs/Bp is a characteristic dimensionless number of the system. We use two plate setups, distinguished by plate length and density. Plate 2 is longer and denser than plate 1. Table 1 details the model parameters used in the simulations. We stress that the models were setup ‘as in the laboratory’ to verify that the simulations had a similar behaviour to that of the laboratory experiments. An additional quantitative comparison of a metrics representative of the long-term evolution and independent of the model dimensions will be shown in Fig. 6b of Section 3. The scaling ensuring similarity between model and nature was as follows: lengths of the model l were scaled so that 1 cm represent 50 km in nature and the timescale was scaled by the time a slab element takes to sink through the model depth |${t^{\text{nature}}\!/ t^{\text{model}} \!=\! (H^{\text{nature}}/H^{\text{model}}) (\Delta \rho_s gd^2\!/ \mu_m )^{\text{model}}/ (\Delta \rho_s gd^2\!/ \mu_m\! )^{\text{nature}}}$|⁠. Hence, 1 s represents 29.5 kyr when the following nature parameters are used: a plate thickness d = 65 km, a plate density contrast Δρs = 100 kg m−3 or Δρs =80 kg m−3, identical to the contrasts used in the simulations, an upper mantle layer H = 400 km and a viscosity μm = 4 × 1020 Pa s. Our simulations typically span a time period of 1500 s, that is 44 Myr. The model velocities scaled to nature are given by |${v^{\text{nature}}=v^{\text{model}} \times (H^{\text{nature}}/H^{\text{model}})( t^{\text{model}}/t^{\text{nature}})}$|⁠. In this study, 1 mm min–1 thus corresponds to 0.28 cm yr–1. Table 1. Numerical simulation parameters. . Variable . Plate setup . Values . Units . Box dimensions  Box height Ly 1, 2 8.8 cm  Box length Lx 1 50 cm 2 56 cm  Box width Lz 1, 2 50 cm Mantle fluid  Gravitational acceleration g 1, 2 9.8 m s−2  Mantle fluid density ρm 1, 2 1417 kg m−3  Mantle fluid viscosity μm 1, 2 110 Pa s  Mantle fluid depth H 1, 2 8 cm Sticky layer  Sticky Layer density ρa 1, 2 1.0 kg m−3  Sticky Layer viscosity μa 1, 2 0.5 Pa s  Sticky Layer thickness h 1, 2 0.8 cm Plate  Plate density ρs 1 1497 kg m−3 2 1517 kg m−3  Plate edge density ρe 1 1427 kg m−3  Plate density contrast |Δρs| = |ρm − ρs| 1 80 kg m−3 2 100 kg m−3  Plate Newtonian viscosity μs 1, 2 57 000 Pa s  Plate thickness d 1, 2 1.3 cm  Plate length (after slab tip initiation) L(0) 1 30 cm 2 36 cm  Plate width W 1, 2 20 cm  Plate edge width |$W_e$| 1, 2 0.45 cm Initial tip  Slab initial tip length l(0) 1, 2 3 cm  Initial dip angle of the slab θ0 1, 2 35 degrees Plume  Plume source location (xp, yp, zp) 1, 2 (12,0,0) cm  Plume radius rp 1, 2 Variable* cm  Plume fluid density ρp 1, 2 1397 kg m−3  Plume fluid viscosity μp 1, 2 5 Pa s  Plume density contrast Δρp = ρm − ρp 1, 2 20 kg m−3  Plume volume flux Q 1, 2 Variable* cm3 s−1 . Variable . Plate setup . Values . Units . Box dimensions  Box height Ly 1, 2 8.8 cm  Box length Lx 1 50 cm 2 56 cm  Box width Lz 1, 2 50 cm Mantle fluid  Gravitational acceleration g 1, 2 9.8 m s−2  Mantle fluid density ρm 1, 2 1417 kg m−3  Mantle fluid viscosity μm 1, 2 110 Pa s  Mantle fluid depth H 1, 2 8 cm Sticky layer  Sticky Layer density ρa 1, 2 1.0 kg m−3  Sticky Layer viscosity μa 1, 2 0.5 Pa s  Sticky Layer thickness h 1, 2 0.8 cm Plate  Plate density ρs 1 1497 kg m−3 2 1517 kg m−3  Plate edge density ρe 1 1427 kg m−3  Plate density contrast |Δρs| = |ρm − ρs| 1 80 kg m−3 2 100 kg m−3  Plate Newtonian viscosity μs 1, 2 57 000 Pa s  Plate thickness d 1, 2 1.3 cm  Plate length (after slab tip initiation) L(0) 1 30 cm 2 36 cm  Plate width W 1, 2 20 cm  Plate edge width |$W_e$| 1, 2 0.45 cm Initial tip  Slab initial tip length l(0) 1, 2 3 cm  Initial dip angle of the slab θ0 1, 2 35 degrees Plume  Plume source location (xp, yp, zp) 1, 2 (12,0,0) cm  Plume radius rp 1, 2 Variable* cm  Plume fluid density ρp 1, 2 1397 kg m−3  Plume fluid viscosity μp 1, 2 5 Pa s  Plume density contrast Δρp = ρm − ρp 1, 2 20 kg m−3  Plume volume flux Q 1, 2 Variable* cm3 s−1 * See Table 2 Open in new tab Table 1. Numerical simulation parameters. . Variable . Plate setup . Values . Units . Box dimensions  Box height Ly 1, 2 8.8 cm  Box length Lx 1 50 cm 2 56 cm  Box width Lz 1, 2 50 cm Mantle fluid  Gravitational acceleration g 1, 2 9.8 m s−2  Mantle fluid density ρm 1, 2 1417 kg m−3  Mantle fluid viscosity μm 1, 2 110 Pa s  Mantle fluid depth H 1, 2 8 cm Sticky layer  Sticky Layer density ρa 1, 2 1.0 kg m−3  Sticky Layer viscosity μa 1, 2 0.5 Pa s  Sticky Layer thickness h 1, 2 0.8 cm Plate  Plate density ρs 1 1497 kg m−3 2 1517 kg m−3  Plate edge density ρe 1 1427 kg m−3  Plate density contrast |Δρs| = |ρm − ρs| 1 80 kg m−3 2 100 kg m−3  Plate Newtonian viscosity μs 1, 2 57 000 Pa s  Plate thickness d 1, 2 1.3 cm  Plate length (after slab tip initiation) L(0) 1 30 cm 2 36 cm  Plate width W 1, 2 20 cm  Plate edge width |$W_e$| 1, 2 0.45 cm Initial tip  Slab initial tip length l(0) 1, 2 3 cm  Initial dip angle of the slab θ0 1, 2 35 degrees Plume  Plume source location (xp, yp, zp) 1, 2 (12,0,0) cm  Plume radius rp 1, 2 Variable* cm  Plume fluid density ρp 1, 2 1397 kg m−3  Plume fluid viscosity μp 1, 2 5 Pa s  Plume density contrast Δρp = ρm − ρp 1, 2 20 kg m−3  Plume volume flux Q 1, 2 Variable* cm3 s−1 . Variable . Plate setup . Values . Units . Box dimensions  Box height Ly 1, 2 8.8 cm  Box length Lx 1 50 cm 2 56 cm  Box width Lz 1, 2 50 cm Mantle fluid  Gravitational acceleration g 1, 2 9.8 m s−2  Mantle fluid density ρm 1, 2 1417 kg m−3  Mantle fluid viscosity μm 1, 2 110 Pa s  Mantle fluid depth H 1, 2 8 cm Sticky layer  Sticky Layer density ρa 1, 2 1.0 kg m−3  Sticky Layer viscosity μa 1, 2 0.5 Pa s  Sticky Layer thickness h 1, 2 0.8 cm Plate  Plate density ρs 1 1497 kg m−3 2 1517 kg m−3  Plate edge density ρe 1 1427 kg m−3  Plate density contrast |Δρs| = |ρm − ρs| 1 80 kg m−3 2 100 kg m−3  Plate Newtonian viscosity μs 1, 2 57 000 Pa s  Plate thickness d 1, 2 1.3 cm  Plate length (after slab tip initiation) L(0) 1 30 cm 2 36 cm  Plate width W 1, 2 20 cm  Plate edge width |$W_e$| 1, 2 0.45 cm Initial tip  Slab initial tip length l(0) 1, 2 3 cm  Initial dip angle of the slab θ0 1, 2 35 degrees Plume  Plume source location (xp, yp, zp) 1, 2 (12,0,0) cm  Plume radius rp 1, 2 Variable* cm  Plume fluid density ρp 1, 2 1397 kg m−3  Plume fluid viscosity μp 1, 2 5 Pa s  Plume density contrast Δρp = ρm − ρp 1, 2 20 kg m−3  Plume volume flux Q 1, 2 Variable* cm3 s−1 * See Table 2 Open in new tab Our models include several simplifications compared to natural systems. First, they are isothermal and isoviscous: the assumed rheological behaviour is independent of both temperature and stress and plumes are compositional. This implies that thermal entrainment into plumes and thermomechanical erosion of the lithosphere or slab by the plume are not considered. Secondly, we constrain the subduction dynamics to only exhibit slab retreat by fixing one end of the plate to the side wall of the domain. In doing so, the plume head propagates beneath the plate without imposed velocity at its top surface as the plate does not have a trench-normal subducting plate velocity (i.e. the subducting velocity is simply the trench-normal trench velocity). Thirdly, our models do not include a low-viscosity asthenospheric layer beneath the plate. Such a low-viscosity layer would enhance the plume head propagation and return flow in the mantle. The simplifications do not impede the interrelationship between the head and conduit dynamics, which we wish to highlight. Snapshots of the system evolution from the time of the plume’s arrival beneath the plate are shown in Appendix A for two representative runs associated with a ‘strong’ and a ‘weak’ plume relative to the plate. The evolution is qualitatively the same as that previously described in laboratory experiments (Mériaux et al. 2015, 2016). 2.2 Numerical method and workflow The dynamics of the slab-plume system described in Section 2.1 is governed by incompressible Stokes flow. The numerical solution of the Stokes flow problem was obtained using the parallel particle-in-cell finite element code Underworld2 (Moresi et al. 2018; Mansour et al. 2020). Underworld2 is the python-friendly version of Underworld (Moresi et al. 2007), which utilises the LavaVu library for inline visualisation (Kaluza et al. 2019). The discretized equations were solved using the multigrid approach with the PETSc package (Balay et al. 1997, 2017a, b). We took advantage of the system symmetry about the hinge-perpendicular z-axis (see Fig. 1) to halve the computational domain in that direction. In all simulations, the numerical grid was uniform in the y direction with a resolution Δy = Ly/Ny = 0.11 cm, where Ny is the number of elements in the y direction. The plate was vertically resolved by 12 elements. In the x- and z-directions, mesh refinement was applied to the numerical grids via mapping functions, resulting in a reduction of the grid cell size in the plate/plume region along x and z. Across all simulations, the minimum and maximum grid resolutions Δx in the x direction ranged from 0.16 to 0.25 cm and 0.40 to 0.61 cm, whereas the minimum and maximum grid resolutions Δz in the z-direction ranged from 0.27 to 0.36 cm and 0.66 to 0.89 cm. Various mesh resolutions were selected upon the requirement that the plume volume Vp(t) that was numerically calculated at each time step and the theoretical plume volume ∫Qdt did not differ by more than  3%. Models with a smaller plume radius and volume flux required meshes with smaller element dimensions, obtained by using more elements in each direction, and/or more aggressive mesh refinement. Typically, the plume source radius rp was resolved by three to six elements in the x and z directions. Two sets of subduction simulations differentiated by a low and a large Bs were run (see Table 2). Table 2. Numerical simulations described in this paper listed by run identification. Slab–plume simulations of each plate setup (2–5 and 7–9) are listed in ascending order of their buoyancy flux ratio. N = Nx × Ny × Nz, where Ni = x, y, z is the number of elements used in the i direction, represents the total number of elements used in the run. Q/2 is the plume volume flux in half the domain. The initial plume head aspect ratio |$\mathcal {H}_0/\mathcal {R}_0$| is determined at the time when the top surface of the plume reaches 6 cm, marking the start of the horizontal spreading. It represents the measure of the plume thickness from the vizualization database at that time divided by |$\mathcal {R}_0=(x_d(t_s)-x_u(t_s))/2$|⁠. Run Id. . Plate setup . N . rp . Q/2 . |$|u_s^s|$| . |$|u_s^r|$| . |$|u_p^r|$| . Bs/Bp . |$\mathcal {H}_0/\mathcal {R}_0$| . . . . (cm) . (cm3 s−1) . (mm min−1) . (mm min−1) . (mm min−1) . . . 1 1 144 × 80 × 96 – – 5.82 6.72 – – – P2 1 144 × 80 × 96 1.1855 0.1520 – – 10.12 0 0.84 2 1 144 × 80 × 96 1.1855 0.1520 5.84 5.65 10.12 3.22 0.96 3 1 160 × 80 × 96 1.0960 0.1111 5.81 6.01 9.09 4.69 0.96 4 1 160 × 80 × 96 0.9216 0.0555 5.83 6.08 7.19 9.49 0.99 P5 1 192 × 80 × 128 0.7870 0.0295 – – 5.91 0 0.87 5 1 192 × 80 × 128 0.7870 0.0295 5.86 5.99 5.91 17.58 0.75 6 2 160 × 80 × 96 – – 7.67 8.79 – – – 7 2 160 × 80 × 96 1.2718 0.2014 7.54 7.39 11.19 3.97 0.85 8 2 160 × 80 × 96 1.1855 0.1520 7.50 7.60 10.05 5.42 0.90 9 2 240 × 80 × 128 0.7870 0.0295 7.66 8.10 5.93 29.72 0.94 Run Id. . Plate setup . N . rp . Q/2 . |$|u_s^s|$| . |$|u_s^r|$| . |$|u_p^r|$| . Bs/Bp . |$\mathcal {H}_0/\mathcal {R}_0$| . . . . (cm) . (cm3 s−1) . (mm min−1) . (mm min−1) . (mm min−1) . . . 1 1 144 × 80 × 96 – – 5.82 6.72 – – – P2 1 144 × 80 × 96 1.1855 0.1520 – – 10.12 0 0.84 2 1 144 × 80 × 96 1.1855 0.1520 5.84 5.65 10.12 3.22 0.96 3 1 160 × 80 × 96 1.0960 0.1111 5.81 6.01 9.09 4.69 0.96 4 1 160 × 80 × 96 0.9216 0.0555 5.83 6.08 7.19 9.49 0.99 P5 1 192 × 80 × 128 0.7870 0.0295 – – 5.91 0 0.87 5 1 192 × 80 × 128 0.7870 0.0295 5.86 5.99 5.91 17.58 0.75 6 2 160 × 80 × 96 – – 7.67 8.79 – – – 7 2 160 × 80 × 96 1.2718 0.2014 7.54 7.39 11.19 3.97 0.85 8 2 160 × 80 × 96 1.1855 0.1520 7.50 7.60 10.05 5.42 0.90 9 2 240 × 80 × 128 0.7870 0.0295 7.66 8.10 5.93 29.72 0.94 Open in new tab Table 2. Numerical simulations described in this paper listed by run identification. Slab–plume simulations of each plate setup (2–5 and 7–9) are listed in ascending order of their buoyancy flux ratio. N = Nx × Ny × Nz, where Ni = x, y, z is the number of elements used in the i direction, represents the total number of elements used in the run. Q/2 is the plume volume flux in half the domain. The initial plume head aspect ratio |$\mathcal {H}_0/\mathcal {R}_0$| is determined at the time when the top surface of the plume reaches 6 cm, marking the start of the horizontal spreading. It represents the measure of the plume thickness from the vizualization database at that time divided by |$\mathcal {R}_0=(x_d(t_s)-x_u(t_s))/2$|⁠. Run Id. . Plate setup . N . rp . Q/2 . |$|u_s^s|$| . |$|u_s^r|$| . |$|u_p^r|$| . Bs/Bp . |$\mathcal {H}_0/\mathcal {R}_0$| . . . . (cm) . (cm3 s−1) . (mm min−1) . (mm min−1) . (mm min−1) . . . 1 1 144 × 80 × 96 – – 5.82 6.72 – – – P2 1 144 × 80 × 96 1.1855 0.1520 – – 10.12 0 0.84 2 1 144 × 80 × 96 1.1855 0.1520 5.84 5.65 10.12 3.22 0.96 3 1 160 × 80 × 96 1.0960 0.1111 5.81 6.01 9.09 4.69 0.96 4 1 160 × 80 × 96 0.9216 0.0555 5.83 6.08 7.19 9.49 0.99 P5 1 192 × 80 × 128 0.7870 0.0295 – – 5.91 0 0.87 5 1 192 × 80 × 128 0.7870 0.0295 5.86 5.99 5.91 17.58 0.75 6 2 160 × 80 × 96 – – 7.67 8.79 – – – 7 2 160 × 80 × 96 1.2718 0.2014 7.54 7.39 11.19 3.97 0.85 8 2 160 × 80 × 96 1.1855 0.1520 7.50 7.60 10.05 5.42 0.90 9 2 240 × 80 × 128 0.7870 0.0295 7.66 8.10 5.93 29.72 0.94 Run Id. . Plate setup . N . rp . Q/2 . |$|u_s^s|$| . |$|u_s^r|$| . |$|u_p^r|$| . Bs/Bp . |$\mathcal {H}_0/\mathcal {R}_0$| . . . . (cm) . (cm3 s−1) . (mm min−1) . (mm min−1) . (mm min−1) . . . 1 1 144 × 80 × 96 – – 5.82 6.72 – – – P2 1 144 × 80 × 96 1.1855 0.1520 – – 10.12 0 0.84 2 1 144 × 80 × 96 1.1855 0.1520 5.84 5.65 10.12 3.22 0.96 3 1 160 × 80 × 96 1.0960 0.1111 5.81 6.01 9.09 4.69 0.96 4 1 160 × 80 × 96 0.9216 0.0555 5.83 6.08 7.19 9.49 0.99 P5 1 192 × 80 × 128 0.7870 0.0295 – – 5.91 0 0.87 5 1 192 × 80 × 128 0.7870 0.0295 5.86 5.99 5.91 17.58 0.75 6 2 160 × 80 × 96 – – 7.67 8.79 – – – 7 2 160 × 80 × 96 1.2718 0.2014 7.54 7.39 11.19 3.97 0.85 8 2 160 × 80 × 96 1.1855 0.1520 7.50 7.60 10.05 5.42 0.90 9 2 240 × 80 × 128 0.7870 0.0295 7.66 8.10 5.93 29.72 0.94 Open in new tab Along the base of the domain, we imposed that the velocity satisfies |$(0,v(r),0)$|⁠. On the left and right boundaries of the domain we prescribed that all components of the velocity were zero (e.g. a no-slip condition), whist on the back boundary we imposed that the normal velocity was zero and the tangential (shear) stress was zero (e.g. a symmetry condition). A thin lubrication layer along the top of the domain was employed to simulate a free surface (Schmeling et al. 2008; Crameri et al. 2012). Along the upper boundary of the domain we imposed zero normal and zero tangential stress (i.e. an open boundary) to accommodate the volume change due to the plume influx. The large viscosity contrasts at the interfaces (plate-air and plate-mantle) were dealt with by using a piecewise constant viscosity within each finite element obtained by harmonically averaging the shear viscosity defined on each material point contained in a given element. Furthermore, following our subduction benchmark study (Mériaux et al. 2018) and the development of plate edge instability, we modified the plate edge by lowering its density ρe over a very small width W compared to the plate half-width |$W/2$| (see Table 1), which in effect retards the instability. Compared to the homogeneous plate, this reduces the initial buoyancy |$\Delta \rho_s W d$| by a factor |$R=1+2W_e/W\left( \Delta \rho_e/\Delta \rho_s -1 \right)$|⁠. In our simulations, R = 0.92 was found to control the plate edge instability without compromising the retreat and sinking rates and plume propagation. For similarity with the laboratory conditions of Mériaux et al. (2015, 2016, 2018), a no-slip velocity condition was also applied to the top of the plate over an inward length of 6 cm from its fixed edge, and over its entire width. Thirty particles carrying the composition (density and viscosity) of the materials were initially distributed in each grid cell, which resulted in a swarm of at least 3.3 × 107 Lagrangian particles. They were advected using a second order Runge-Kutta scheme. In addition, passive tracers were initially randomly distributed at the height of the plate base. Other passive tracers were initially aligned within the plate interior in a series of 10 along-plate lines equally spaced across the plate half-width (see Fig. 1). Those tracers were used to image the slab-induced flow and estimate the plate hinge retreat xh(t) as part of the simulation time-stepping workflow. Also, as part of the simulation time-stepping workflow, the height of the slab tip yt(t), the plume head horizontal maximum position beneath the plate moving toward the trench xu(t), the plume head horizontal maximum position beneath the plate moving away from the trench xd(t), and the plume head maximum position beneath the plate and parallel to the trench zl(t) were calculated from the Lagrangian swarm. The time-stepping workflow also included a calculation of the horizontal position of the plume centre of mass xc(yi, t) in the plane z = 0 within 26 horizontal layers evenly distributed between the base of the box and the base of the plate. For some simulations (7–9), the maximum extent of the plume material within each layer in the z-direction |$z_{max}(y_i,t)$| was also estimated. All these positions are shown in Fig. 1. Although the last two of these markers give some information on the plume head, they are primarily used to track the plume conduit evolution as shown in Fig. 2. Ultimately, we identify the positions xc(yi, t) and |$z_{max}(y_i,t)$| at the specific height yi marking the final height of the conduit before the horizontal plume spreading by finding where the derivatives of xc(yi, t) and |$z_{max}(y_i,t)$| with height most significantly change. In the rest of this paper, for simplicity, we will refer to these ‘feeding’ positions as being xc(t) and |$r_c^z$|⁠, a shown in Fig. 2. Figure 2. Open in new tabDownload slide (a) Position of the plume centre of mass xc(yi, t) in the x direction and plane z = 0 shown at several times t in run 9. (b) Maximum extent of the plume in the z direction |$z_{max}(y_i,t)$| for each layer i of mid-height yi shown at several times t in run 9. The horizontal black dashed lines indicate the base of the plate y = H − d. The filled circles mark the vertical end of the conduit. In (a), xc is the position of the conduit centre of mass at the base of the plume head. In (b) |$r_c^z$| is the conduit extent measured at the base of the plume head. Figure 2. Open in new tabDownload slide (a) Position of the plume centre of mass xc(yi, t) in the x direction and plane z = 0 shown at several times t in run 9. (b) Maximum extent of the plume in the z direction |$z_{max}(y_i,t)$| for each layer i of mid-height yi shown at several times t in run 9. The horizontal black dashed lines indicate the base of the plate y = H − d. The filled circles mark the vertical end of the conduit. In (a), xc is the position of the conduit centre of mass at the base of the plume head. In (b) |$r_c^z$| is the conduit extent measured at the base of the plume head. In our configuration, the region of interest is at the rear of the slab where the plume source has been located. Fig. 3 shows the typical flow field in the presence of a plume at the rear of a slab. The presence of the plume limits the subduction-induced poloidal flow to the region between the slab and the plume, while the toroidal flow caused by subduction opposes the plume head buoyancy-driven spreading close to the plane of symmetry z = 0 but the two flows merge at a distance from this plane of about one third of plane width and, at a position in the x direction either ahead of the plume head upstream front xu < x < xh (‘strong’ plume) or as far back as the plume head centre xd < x < xu (‘weak’ plume) (see Figs 3a and 3b, respectively). With respect to the velocity components, those of specific interest in the region of slab–plume interaction beneath the plate are |$v_x$| and |$v_z$|. Our post-processing workflow will thus include the extraction of further quantitative information on |$v_x$| and |$v_z$| through velocity profiles made with and without plume. Figure 3. Open in new tabDownload slide Side and top views of the flow field in (a) run 2 at t = 598 s (Bs/Bp = 3.22) and (b) run 5 at t = 807 s (Bs/Bp = 17.58). These two times correspond to the same slab-plume distance xh − xu = 7 cm. The views are in the planes z = 0 cm and y = 5.3 cm, respectively. y = 5.3 cm is located at approximately mid-height in the plume head. Figure 3. Open in new tabDownload slide Side and top views of the flow field in (a) run 2 at t = 598 s (Bs/Bp = 3.22) and (b) run 5 at t = 807 s (Bs/Bp = 17.58). These two times correspond to the same slab-plume distance xh − xu = 7 cm. The views are in the planes z = 0 cm and y = 5.3 cm, respectively. y = 5.3 cm is located at approximately mid-height in the plume head. 3 INFLUENCE OF THE PLUME ON SUBDUCTION The influence of the plume on the subduction dynamics is assessed through the analysis of the slab’s sinking and retreat without and with a plume. Figs 4 and 5 show the vertical position yt of the slab tip at the slab centreline and the horizontal position of the plate hinge to the plate fixing wall xh as a function of time t for all the simulations. The slab’s sinking is marked by three successive phases characteristic of the initial slab bending (phase I), the sinking phase which includes the steady sinking followed by a slow down (phase II) and its resting phase on the bottom surface (phase III). The same phases I–III related to bending and pre- and post-sinking can be identified in the retreat. However, in the presence of a plume, the retreat further slows down once the front of the plume head xu reaches the slab (phase IV). All plumes influence the retreat across the different phases, except the bending phase I, which is quasi free from plume influence. Figure 4. Open in new tabDownload slide Vertical position of the slab tip yt measured at the slab centreline as a function of time t. The inset shows the four phases I–IV typically identified throughout the slab retreat. The solid and dashed lines distinguish runs by initial plate setups, respectively, 1 and 2 (see Table 2). Plate 2 is longer and denser than plate 1 (see Table 1). Figure 4. Open in new tabDownload slide Vertical position of the slab tip yt measured at the slab centreline as a function of time t. The inset shows the four phases I–IV typically identified throughout the slab retreat. The solid and dashed lines distinguish runs by initial plate setups, respectively, 1 and 2 (see Table 2). Plate 2 is longer and denser than plate 1 (see Table 1). Figure 5. Open in new tabDownload slide Horizontal position of the hinge of the plate xh to the plate fixing wall measured at the slab centreline as a function of time t. The inset shows the four phases I–IV typically identified throughout the slab retreat. The solid and dashed lines distinguish runs by initial plate setups, respectively, 1 and 2 (see Table 2). Plate 2 is longer and denser than plate 1 (see Table 1). Figure 5. Open in new tabDownload slide Horizontal position of the hinge of the plate xh to the plate fixing wall measured at the slab centreline as a function of time t. The inset shows the four phases I–IV typically identified throughout the slab retreat. The solid and dashed lines distinguish runs by initial plate setups, respectively, 1 and 2 (see Table 2). Plate 2 is longer and denser than plate 1 (see Table 1). To quantify the plume influence on subduction, we first estimated the sinking rates |$|u_s^s|$| by linear fits of the slab tip vertical position yt(t) over a depth interval in the steady sinking phase between y = 3.5 cm and y = 1.5 cm for all simulations (Table 2). For the retreat, we evaluated a mean retreat rate |$|u_s^r|$| by taking the average of the retreat velocities calculated as the first-order derivative of the hinge horizontal position xh(t) with respect to time over phase II–IV (included). Such values of |$|u_s^r|$| are representative of the overall retreat rate, which is used to scale the slab buoyancy ratio (Table 2) and assess the general influence of the plume on the slab retreat (Fig. 6a). We further evaluated the retreat rates in phase IV relative to phase III to get an insight into the ultimate influence of the plume once at the slab (Fig. 6b). These last estimates were compared to previous similar laboratory experiments (Mériaux et al. 2015). Figure 6. Open in new tabDownload slide Influence of mantle plumes on retreat rate as a function of buoyancy flux ratio Bs/Bp. (a) Change of retreat rate of the slab–plume simulations |$(u_s^r) _{SP}$| relative to slab-only simulation |$(u_s^r )_{SO}$| for each plate setting, |$|\Delta (u_s^r)_{SP/SO}|= |((u_s^r) _{SP} - (u_s^r)_{SO}) / (u_s^r )_{SO}|$|⁠. (b) Change of retreat rate of the slab–plume simulations in the final phase IV relative to phase III after the slab has reached the bottom. The change is defined as |$|\Delta u_s^r|_{IV/III}=|((u_s^r) _{III} - (u_s^r)_{IV}) / (u_s^r )_{III}|$| (blue diamond). The laboratory data (red circle) are from Mériaux et al. (2015). The shaded areas corresponds to simulations and experiments, which have seen plume material escape from the plate at the free edge. Figure 6. Open in new tabDownload slide Influence of mantle plumes on retreat rate as a function of buoyancy flux ratio Bs/Bp. (a) Change of retreat rate of the slab–plume simulations |$(u_s^r) _{SP}$| relative to slab-only simulation |$(u_s^r )_{SO}$| for each plate setting, |$|\Delta (u_s^r)_{SP/SO}|= |((u_s^r) _{SP} - (u_s^r)_{SO}) / (u_s^r )_{SO}|$|⁠. (b) Change of retreat rate of the slab–plume simulations in the final phase IV relative to phase III after the slab has reached the bottom. The change is defined as |$|\Delta u_s^r|_{IV/III}=|((u_s^r) _{III} - (u_s^r)_{IV}) / (u_s^r )_{III}|$| (blue diamond). The laboratory data (red circle) are from Mériaux et al. (2015). The shaded areas corresponds to simulations and experiments, which have seen plume material escape from the plate at the free edge. Compared to run 1, the increase of the plate density contrast |Δρs| by 20 % in run 6 results in an increase of |$|u_s^s|$| and |$|u_s^r|$| by 30 % (see Table 2). Consequently, slab sinking and retreat rates are higher for plate setup 2 than for plate setup 1. For each given plate setup, the sinking rates |$|u_s^s|$| are similar between simulations (see Fig. 4 and Table 2). By contrast, retreat rates |$|u_s^r|$| vary as a function of plume flux Q as shown in Fig. 5 and Table 2. Overall, the influence of the plume on the retreat rates increases with Bs/Bp and varies from 7 % to 15 % between our plume–slab simulations as shown in Fig. 6a. When the slab–plume retreat starts its phase III, which differentiates itself from its slab-only counterpart due to the remote plume influence, the slab’s hinge distance from the plume source xh − xp ranges from 8 to 17 cm. We note that the slab-plume distance at the onset of phase III in run 7 was 14.5 cm, which is 36 % longer than the slab–plume distance at the onset of phase III in run 2 with a similar buoyancy (9.2 cm). The ultimate direct impact of the plume on the slab further affects the slab retreat, whose trend is to decrease as Bs/Bp increases as shown in Fig. 6b. At small ratios Bs/Bp, there is a variability in the slab response to the plume impact, which seems to coincide with the escape of plume material from beneath the plate at the free edge during phase IV. In runs 7 and 8, this leakage of plume buoyancy beneath the plate, when the extent of plume spreading in the z-direction exceeded the plate width, occurred at times t = 1217 and 1379 s, respectively. The simulations and laboratory data are consistent with one another. 4 INTER-RELATIONSHIP BETWEEN THE PLUME HEAD AND ITS CONDUIT Two processes characterize the plume dynamics: the deflection and deformation of the conduit and the non-axisymmetric plume head spreading under the plate and ultimately beneath the slab. The horizontal deflection of the conduit and its enlargement in the direction parallel to the trench as a function of height can be seen in Fig. 2. The plume head asymmetric spreading shown in Fig. 1 will be further detailed in Fig. 16. In this section, the aim is to show that these two processes are interrelated. To show this interrelationship, we use as a reference the ambient flow generated by subduction only. 4.1 Mantle flow at the rear of the slab We first discuss the ambient mantle flow in the absence of a plume (Figs 7 a and 8a), the impact of the plume on the mantle flow will be highlighted in Sections 4.2 and 4.3. As shown in Fig. 7a for a slab-only simulation, a narrow weak and almost constant shearing layer in the x direction is located beneath the plate between y = 6.4 cm and y = 6.7 cm. The plume head beneath the plate is thus to be slightly sheared at its upper surface towards the slab. Below, the shear velocity depth profile |$v_x(y)$| in the x direction increases. It is non-steady through time, but it remains well represented by a parabolic function of depth in the interval y = [0, 6.4] cm with its peak at y = 3.2 cm. This shear is to oppose the upstream head front spreading and deflect the plume conduit away from the slab. Further downstream, the magnitude of |$v_x$|decreases (profile not shown). For instance, in the plane z = 0, |$v_x$| decreased by a factor ∼6 between x =18 cm and 6 cm at the level of the head (y = 5 cm). Shearing velocity component |$v_z(y)$| in the z direction has a similar parabolic depth profile but it is about nine times smaller in magnitude, that is |$|v_x(y)| \sim 9 |v_z(y)|$|⁠. However, this trend reverses below the plate as the shear in the z direction increases with distance from the plane of symmetry z = 0 (see Fig. 8a). At the plate edge z = 10 cm, |$|v_z(z)| \sim 4 |v_x(z)|$|⁠. This shear velocity component is to contribute to shearing the plume head away from the plane of symmetry z = 0. Figure 7. Open in new tabDownload slide Profiles of the velocity component |$v_x$| (solid curve) and |$v_z$| (dashed curves) at a fixed position x = xp + 6 = 18 cm and z = 1 cm, 4–7 cm at the rear of the slab, as function of height y for three positions of the trench xh =22, 23 and 24 cm. The black horizontal line marks the height of the base of the plate, y = H − d. (a) Slab-only run 1. (b) Slab–plume run 2. Note that in (b) the maximum of |$|v_x|$| is increased by 17 % compared to a) when xh ∼22 cm. On the right, the two sets of 3 panels each correspond from top to bottom to the three positions of the trench xh = 24, 23 and 22 cm. The green line indicates where the profile was made. Figure 7. Open in new tabDownload slide Profiles of the velocity component |$v_x$| (solid curve) and |$v_z$| (dashed curves) at a fixed position x = xp + 6 = 18 cm and z = 1 cm, 4–7 cm at the rear of the slab, as function of height y for three positions of the trench xh =22, 23 and 24 cm. The black horizontal line marks the height of the base of the plate, y = H − d. (a) Slab-only run 1. (b) Slab–plume run 2. Note that in (b) the maximum of |$|v_x|$| is increased by 17 % compared to a) when xh ∼22 cm. On the right, the two sets of 3 panels each correspond from top to bottom to the three positions of the trench xh = 24, 23 and 22 cm. The green line indicates where the profile was made. Figure 8. Open in new tabDownload slide Profiles of the velocity component |$v_x$| (solid curve) and |$v_z$| (dashed curves) at a fixed position x = xp + 6 = 18 cm and y = 6 cm, 4–7 cm at the rear of the slab, as function of z for three positions of the trench xh = 22, 23 and 24 cm. The black horizontal line marks the limit of the plate width. (a) Slab-only run 1. (b) Slab–plume run 2. On the right, the two sets of 3 panels each correspond from top to bottom to the three positions of the trench xh = 24, 23 and 22 cm. The green line indicates where the profile was made. Figure 8. Open in new tabDownload slide Profiles of the velocity component |$v_x$| (solid curve) and |$v_z$| (dashed curves) at a fixed position x = xp + 6 = 18 cm and y = 6 cm, 4–7 cm at the rear of the slab, as function of z for three positions of the trench xh = 22, 23 and 24 cm. The black horizontal line marks the limit of the plate width. (a) Slab-only run 1. (b) Slab–plume run 2. On the right, the two sets of 3 panels each correspond from top to bottom to the three positions of the trench xh = 24, 23 and 22 cm. The green line indicates where the profile was made. 4.2 Plume heads impact on the conduit The conduit deflection is first linked to both the slab- and plume-induced velocity in the mantle at the rear of the slab. Due to incompressibility, the plume head spreading adds to the slab-only return poloidal flow, which is a parabolic function of height, increasing the horizontal shear velocity and lowering the depth of its peak. For instance, in the early stages of the plume head propagation in run 2 and at the upstream position x = xp + 6 cm, the magnitude of the horizontal shear velocity |$v_x(y)$| in the mantle increases by 17 % and its peak is at y =2.7 cm (i.e. a change of 0.5 cm or 25 km compared to the slab-only simulation) as shown in Fig. 7a. The plume head thus contributes to shearing its own conduit. As shown in Fig. 9, this occurs regardless of the value of the buoyancy flux ratio Bs/Bp. Figure 9. Open in new tabDownload slide Increase of the maximum shearing component |$|v_x^{max}|$| in each slab-plume run (SP in y label) relative to its slab-only equivalent run (SO in y label). The measurement was performed systematically through profiles in the plane z = 0 and positions x located 4 cm behind trench positions xh that were determined by the slab-plume runs when xh − xu = 4 cm (i.e. the upstream plume head front is 4 cm away from the slab). The increase varies between 15 and  20%. Figure 9. Open in new tabDownload slide Increase of the maximum shearing component |$|v_x^{max}|$| in each slab-plume run (SP in y label) relative to its slab-only equivalent run (SO in y label). The measurement was performed systematically through profiles in the plane z = 0 and positions x located 4 cm behind trench positions xh that were determined by the slab-plume runs when xh − xu = 4 cm (i.e. the upstream plume head front is 4 cm away from the slab). The increase varies between 15 and  20%. As demonstrated in Appendix B (eqs B1–B3), in steady-state, a parabolic ambient shear profile |$v(y)$| of the form shown in Fig. 7a leads to a shape of the plume conduit that is a cubic function of height over the mantle depth. At a late stage, when the plume head upstream front xu is close to the slab, the slab–plume shear profile departs from a parabolic profile over the heights y = [0, H − d], where H − d is the maximum height of mantle below the plate. This is due to the propagation of the plume head in the interval yb < y < H − d, where yb is at the base of the plume head (see Fig. 10a). As a result, the slab–plume shear profile is better fitted by a cubic function of height of the form given by eq. (B1), where the depth of mantle fluid is yb (see Fig. 10a). The conduit trajectory calculated from the instantaneous slab-only- and slab–plume induced shear profile using eq. (B2) are shown in Fig. 10b. The difference between the computed reference conduit positions and the reconstructed positions was evaluated via the mean absolute errors MAEs as follows: |$MAE=\sum _{i=1}^{N} | x_i^{n}-x_i^{p}| /N$|⁠, where N is the number of points, |$x_i^{n}$| is the numerical horizontal conduit position at height 0 ≤ yi < yb cm and |$x_i^{p}$| is the horizontal conduit position at the same height yi calculated from either the slab-only induced velocity profile or combined slab–plume profile. The MAEs between the computed reference conduit path and the two reconstructed conduit paths from the instantaneous slab-only and slab–plume shear profiles are 0.19 and 0.09, respectively. The slab–plume induced shear profile leads to a better fit of the conduit trajectory but the conduit trajectory calculated from the slab-only instantaneous velocity profile can be used as an indicator of the tilt of the conduit. The fact that calculations from the instantaneous velocities appear reasonable in run 2 (Bs/Bp = 3.22) and runs at low buoyancy flux ratio (Bs/Bp ≤ 9.5) shows that the plume’s response time τ is sufficiently short compared to the timescale associated with the evolution of the ambient shear. We note that τ is expected to vary as the inverse of the plume Stokes velocity |$v_S$| and increase with height as |$\tau = y/v_S \propto y/r_p^2$| (Richards & Griffiths 1988). In run 2, τ < 332s. For weaker plumes (e.g., runs 5 and 9; Bs/Bp > 9.5), the agreement between the reconstructed and simulated conduit paths is not as good in the last centimetres of the mantle layer (⁠|$y\gt 30-40\, {\rm \%}\, (H-d)=$| 4–4.70 cm), where the conduit response time to the ambient is expected to be longest due to both the heights y and the smaller values of rp. For instance, the only effect of rp implies that |$\tau_{\text{run 5}}\,=\, 2.3 \tau_{\text{run 2}}$|⁠. This difference in plume deflection is further outlined in Fig. 11, where the efficiency of the estimate for the conduit offset derived from eqs (B1) and (B4) is evaluated as a function of buoyancy number Bs/Bp. The simulated conduit offset x(y) − xp is measured at a height of y = (H − d)/5 = 4 cm (i.e. close to the plume head base for all simulations) and compared to the slab-only reconstructed offset given by $$\begin{eqnarray*} \left|x \left(\frac{3(H-d)}{5}\right)-x_p\right|=\frac{27(H-d)}{250} \frac{4v_x^{max}}{v_S}, \end{eqnarray*}$$(6) where |$v_S$| is the plume Stokes velocity as defined in Appendix B and |$v_x^{max}$| is the maximum horizontal shear of the slab-only simulation for a trench position xh, which matches that of the slab–plume simulation when xh − xu = 2 cm. All other symbols are defined in Table 1 and Fig. 1. The deflection analytically calculated from the instantaneous slab-only profile gives a reasonable estimate when Bs/Bp ≤ 9.5. Figure 10. Open in new tabDownload slide Analysis of the conduit trajectory when the plume head upstream front xu is close to the slab (xh − xu = 1.2 cm; xh = 19 cm) in run 2. (a) Shear velocity depth profiles and their polynomial fits at x = 15 cm and z = 0 cm in slab-only run 1 and slab–plume run 2 when xh = 19 cm. The horizontal dashed line marks the base of the propagating plume head. (b) Numerical conduit trajectory, fitted trajectory and predicted trajectories from slab-only and slab–plume shear velocity profile when xh = 19 cm. Figure 10. Open in new tabDownload slide Analysis of the conduit trajectory when the plume head upstream front xu is close to the slab (xh − xu = 1.2 cm; xh = 19 cm) in run 2. (a) Shear velocity depth profiles and their polynomial fits at x = 15 cm and z = 0 cm in slab-only run 1 and slab–plume run 2 when xh = 19 cm. The horizontal dashed line marks the base of the propagating plume head. (b) Numerical conduit trajectory, fitted trajectory and predicted trajectories from slab-only and slab–plume shear velocity profile when xh = 19 cm. Figure 11. Open in new tabDownload slide Deflection of the conduit relative to the plume source |x(y) − xp| as a function of buoyancy flux ratio Bs/Bp. The deflection is measured at height y = 3(H − d)/5 = 4 cm, which corresponds to three-fifths of the mantle layer thickness below the plate. This specific height has been chosen for reference for all simulations as it is nearby the base of all plumes. It is also systematically measured at a time when the plume head upstream front xu is near the slab (xh − xu = 2 cm) and for a matching xh position in the corresponding slab-only runs 1 & 6. The deflection measured in the simulations is compared with the analytic solution eq. (6) derived from eqs (B1) and (B4). The maximum slab-only ambient shear velocity |$v_x^{max}$| has been determined from shear depth profiles made at x = 15 cm and z = 0 (i.e. in the same way as shown in Fig. 10). Figure 11. Open in new tabDownload slide Deflection of the conduit relative to the plume source |x(y) − xp| as a function of buoyancy flux ratio Bs/Bp. The deflection is measured at height y = 3(H − d)/5 = 4 cm, which corresponds to three-fifths of the mantle layer thickness below the plate. This specific height has been chosen for reference for all simulations as it is nearby the base of all plumes. It is also systematically measured at a time when the plume head upstream front xu is near the slab (xh − xu = 2 cm) and for a matching xh position in the corresponding slab-only runs 1 & 6. The deflection measured in the simulations is compared with the analytic solution eq. (6) derived from eqs (B1) and (B4). The maximum slab-only ambient shear velocity |$v_x^{max}$| has been determined from shear depth profiles made at x = 15 cm and z = 0 (i.e. in the same way as shown in Fig. 10). In addition to tilting, the conduit gradually enlarges with height and time in the direction parallel to the trench (z-direction) as shown in Fig. 2b. The arrival of the plume head beneath the slab is synchronous with further enlargement as shown in Fig. 12. Such flaring of the conduit depends on the buoyancy flux ratio. At low buoyancy flux ratios (i.e. relatively strong plumes) like in run 7 or 8, the flow in the conduit is primarily vertical as tilting is limited, whereas the mantle toroidal slab- and plume-induced crossflow stretch the conduit in the z-direction (Fig. 12b). Similar to the deformation of a viscous blob in pure shear flow first described by Taylor (1934), the stretching is exponential (see further details in the Appendix C). The conduit extent at the base of the plume head |$r_c^z(t)$| (see Fig. 2) varies exponentially with time and can be fitted with a function of the form given by eq. (C2) in which the strain rate is on the order of the gradient of |$v_z$| over the plate width at the height of the plume spreading (see Fig. 8). For instance, in run 7, |$r_c^z(t)$| fits the functional |$r_c^z(t)/r_p=1.245 \, \mbox{e}^{0.0003625t}$|⁠, implying that the strain rate |$\dot{\varepsilon }$| was 0.0003625 s−1, that is 0.0217 min−1. Such a value of the strain rate is consistent with that of the gradient of |$v_z$| beneath the plate given by |$\dot{\varepsilon }\sim dv_z/dz \sim v_z/W \sim 0.33/10 =0.033$| min−1, where the value |$v_z=3.3$| mm min–1 is determined from a profile of |$v_z$| similar to that shown in Fig. 8a but for the slab-only run 6. At large buoyancy flux ratio (i.e. relatively weak plumes - cf. run 9) the conduit extent parallel to trench at the base of the plume head |$r_c^z(t)$| does not seem to solely deform viscously as the fitted exponential function, |$r_c^z(t)/r_c=0.4431 \, \mbox{e}^{0.0008671t}$| does not satisfy the limit |$r_c^z(t_0)/r_p \sim 1$|⁠. Furthermore, the value of the strain rate is not consistent with the gradient of |$v_z$| over the plate width at the height of plume spreading. This is likely due to the fact that conduits of weak plumes are more deflected and consequently the laminar flow in the conduit is partitioned into two components: a diminished vertical component and a new horizontal component causing the fluid to flow in the same direction as the mantle flow crossflow. The increase in |$r_c^z$| is accompanied by a decrease of the conduit extent |$r_c^x$| in the x direction to accommodate the constant source flux. |$r_c^x$| is calculated using eq. (C7) assuming that the deformed conduit is elliptic, which is stricto sensu an approximation as the mantle flow on either side of the conduit is not symmetric. By the time the plume head has reached the slab, |$r_c^z(t)$| has increased by 50 % whereas |$r_c^x(t)$| has decreased by  25%. Figure 12. Open in new tabDownload slide (a) Distance between the plume head front xu(t) and the rear of the plate xh(t) − d as function of time t in runs 7–9. When xh(t) − d − xu(t) = 0, the plume material has reached the slab. (b) Length of the conduit |$r_c^z(t)$| (coloured solid curves) in the z direction at the base of the plume head as a function of time t in runs 7–9. |$r_c^z(t)$| has been made dimensionless with the initial conduit radius rp. The dashed black lines are exponential fits. The coloured dashed lines are the length of the conduit |$r_c^x(t)$| in the x direction at the base of the plume head calculated using eq. (C7) and also made dimensionless with rp in runs 7–9. (c) Distance between the front of the current zl(t) as function of time t. In (b) and (c), the diamonds indicate the times at which xh(t) − d − xu(t) = 0 (i.e. when the plume head has reached the slab). Note that the different initial times correspond to the times of plume arrival beneath the surface and start of horizontal spreading. Figure 12. Open in new tabDownload slide (a) Distance between the plume head front xu(t) and the rear of the plate xh(t) − d as function of time t in runs 7–9. When xh(t) − d − xu(t) = 0, the plume material has reached the slab. (b) Length of the conduit |$r_c^z(t)$| (coloured solid curves) in the z direction at the base of the plume head as a function of time t in runs 7–9. |$r_c^z(t)$| has been made dimensionless with the initial conduit radius rp. The dashed black lines are exponential fits. The coloured dashed lines are the length of the conduit |$r_c^x(t)$| in the x direction at the base of the plume head calculated using eq. (C7) and also made dimensionless with rp in runs 7–9. (c) Distance between the front of the current zl(t) as function of time t. In (b) and (c), the diamonds indicate the times at which xh(t) − d − xu(t) = 0 (i.e. when the plume head has reached the slab). Note that the different initial times correspond to the times of plume arrival beneath the surface and start of horizontal spreading. 4.3 Plume feeder conduits impact on plume head spreading The main characteristics of the plume head spreading once the plume has reached the plate are outlined in Fig. 13. Upstream, the propagation front reaches a stagnation point. The time at which this occurs shortens as the slab/plume buoyancy ratio increases. By contrast, the propagation downstream and in the direction parallel to the trench are not impeded. The end-result of the propagation of the fronts can be reduced to the aspect ratio of the head as shown in Fig. 13c. At the time the slab is close to the head, the head’s extent parallel to the trench |$z_l(t)$| is ∼ 1.5 times larger than its length perpendicular to the trench |$w(t)=x_u(t)-x_d(t)$|⁠. Figure 13. Open in new tabDownload slide Plume head spreading. (a) and (b) Propagation of the plume head fronts xu(t), xd(t) andzl(t) as a function of time t in runs 2 (Bs/Bp = 3.22) and 9 (Bs/Bp = 29.72), respectively. The horizontal black dashed lines intercept the stagnation point. (c) Aspect ratio of the plume head measured as zl(t)/(w(t)/2), where w(t) = (xu(t) − xd(t)) in all simulations. Figure 13. Open in new tabDownload slide Plume head spreading. (a) and (b) Propagation of the plume head fronts xu(t), xd(t) andzl(t) as a function of time t in runs 2 (Bs/Bp = 3.22) and 9 (Bs/Bp = 29.72), respectively. The horizontal black dashed lines intercept the stagnation point. (c) Aspect ratio of the plume head measured as zl(t)/(w(t)/2), where w(t) = (xu(t) − xd(t)) in all simulations. The first explanation for the non-axisymmetric spreading of the head is the external mantle flow at the rear of the slab and level of the head (see Section 4.1). In particular, velocities of the plume and the mantle upward in the x direction are opposed as shown in Figs 7b and 8b by the change of direction of |$v_x$| in the presence of a plume. Eventually, it leads to a halt of the plume front and mantle in this direction when the two velocities are equal. The momentum is then transferred in the z-direction. Velocities of the mantle in the z-direction increase in the presence of the plume by a factor ∼2 (see Figs 7b and 8b). This is to contribute driving the plume head in the direction parallel to the trench while decreasing upstream head spreading. Yet, such rationale for the head asymmetric spreading is drawn in complete ignorance of the feeding conditions. When considering the feeding conditions, the first observation is the motion of the feeding source downstream (i.e. in the same direction as the retreat), which accompanies the conduit tilting. For instance, in run 9 (Bs/Bp ∼ 30), the conduit source moved in the x direction away from the slab by a maximum of 5.42 cm (271 km) in 1091 s (32 Myr). The motion of the feeding conduit decreases when the plume buoyancy increases but it is never fully hindered as shown in Fig. 14. Such motion is another factor adding to the asymmetric spreading. Figure 14. Open in new tabDownload slide Evolution of the position of the conduit centre of mass at the base of the plume head xc(t), as shown in Fig. 2, as a function of time t in all slab–plume simulations. Figure 14. Open in new tabDownload slide Evolution of the position of the conduit centre of mass at the base of the plume head xc(t), as shown in Fig. 2, as a function of time t in all slab–plume simulations. Furthermore, at the release by the conduit of the plume volume flux Q upon its arrival beneath the plate, the axial flow in the conduit turns into a radial horizontal outflow into the head. In the absence of external flow, the axial flow in the conduit is entirely vertical and consequently the radial redistribution of the plume flux into the head is equal among all azimuths (Q/2π). Here, however, the tilt and flare of the conduit modify the axial flow in the conduit by generating a horizontal velocity component in the directions of tilt and flare. As a result, when released, the plume flux is no longer equally redistributed into the head over all azimuths; it preferentially goes in the directions of tilt and flare. Importantly, this redistribution is time-dependent as it follows the evolution of the conduit, which means that the source flux condition in the head takes the form |$Q\mathcal {F}(t,\theta )$|⁠. To obtain insight into the latter, we compared in runs 2 and 5 the plume flux released from the conduit in the head upstream (i.e. opposite to the direction of retreat), downstream in the direction of the conduit tilt and retreat and in the z-direction of flaring (see Fig. 2). Fig. 15 shows the positioning of the control volume used to estimated the plume volume flux in each of these three directions. As detailed in Table 3, the plume flux released into the head goes with time preferentially in the z-direction and then downstream compared with upstream. This is even more pronounced in the case of a weak plume. Figure 15. Open in new tabDownload slide Positions of the control volumes used to calculate the volume fluxes released into the head in three directions: upstream, downstream and parallel to the trench (flaring direction; seen from the bottom top view only). Figure 15. Open in new tabDownload slide Positions of the control volumes used to calculate the volume fluxes released into the head in three directions: upstream, downstream and parallel to the trench (flaring direction; seen from the bottom top view only). Table 3. Volume fluxes in the plume head |$q_i=\int v_i dv/ \Delta x_i$|⁠, where i refers to the x or z direction, for runs 2 and 5 at an ‘early’ and evolved stage. The volume fluxes in the head at the earlier time t = 478 s in run 2 are shown as the plume shows little asymmetry. Two control volumes are diametrically opposed in the x-direction while the third is along z. The positions of the control volumes are shown in Fig. 15. The numerical calculation of the control volume |$V_{control}=\int dv$| was within 1.2  % of the prescribed control volume. Run Id. . Bs/Bp . t . Vcontrol . . . . . |$|q_x^{down}| /(Q/2)$| . |$q_z^{lat}/(Q/2)$| . . . (s) . (cm3) . (cm3s−1) . (cm3s−1) . (cm3s−1) . ( %) . ( %) . ( %) . 2 3.65 478 4.4 0.0305 −0.0439 0.0443 20 29 29 2 3.65 959 8.8 0.0217 −0.0374 0.0381 14 24 25 5 20.04 1198 1.8 0.0025 −0.0077 0.0084 8 26 28 5 20.04 1372 1.4 0.00151 −0.0067 0.00792 5 23 27 Run Id. . Bs/Bp . t . Vcontrol . . . . . |$|q_x^{down}| /(Q/2)$| . |$q_z^{lat}/(Q/2)$| . . . (s) . (cm3) . (cm3s−1) . (cm3s−1) . (cm3s−1) . ( %) . ( %) . ( %) . 2 3.65 478 4.4 0.0305 −0.0439 0.0443 20 29 29 2 3.65 959 8.8 0.0217 −0.0374 0.0381 14 24 25 5 20.04 1198 1.8 0.0025 −0.0077 0.0084 8 26 28 5 20.04 1372 1.4 0.00151 −0.0067 0.00792 5 23 27 Open in new tab Table 3. Volume fluxes in the plume head |$q_i=\int v_i dv/ \Delta x_i$|⁠, where i refers to the x or z direction, for runs 2 and 5 at an ‘early’ and evolved stage. The volume fluxes in the head at the earlier time t = 478 s in run 2 are shown as the plume shows little asymmetry. Two control volumes are diametrically opposed in the x-direction while the third is along z. The positions of the control volumes are shown in Fig. 15. The numerical calculation of the control volume |$V_{control}=\int dv$| was within 1.2  % of the prescribed control volume. Run Id. . Bs/Bp . t . Vcontrol . . . . . |$|q_x^{down}| /(Q/2)$| . |$q_z^{lat}/(Q/2)$| . . . (s) . (cm3) . (cm3s−1) . (cm3s−1) . (cm3s−1) . ( %) . ( %) . ( %) . 2 3.65 478 4.4 0.0305 −0.0439 0.0443 20 29 29 2 3.65 959 8.8 0.0217 −0.0374 0.0381 14 24 25 5 20.04 1198 1.8 0.0025 −0.0077 0.0084 8 26 28 5 20.04 1372 1.4 0.00151 −0.0067 0.00792 5 23 27 Run Id. . Bs/Bp . t . Vcontrol . . . . . |$|q_x^{down}| /(Q/2)$| . |$q_z^{lat}/(Q/2)$| . . . (s) . (cm3) . (cm3s−1) . (cm3s−1) . (cm3s−1) . ( %) . ( %) . ( %) . 2 3.65 478 4.4 0.0305 −0.0439 0.0443 20 29 29 2 3.65 959 8.8 0.0217 −0.0374 0.0381 14 24 25 5 20.04 1198 1.8 0.0025 −0.0077 0.0084 8 26 28 5 20.04 1372 1.4 0.00151 −0.0067 0.00792 5 23 27 Open in new tab The asymmetry in the propagation of the plume head relative to the source as shown in Fig. 16 is thus the result of the combined feeding conditions and ambient flow in the mantle. These spreading conditions greatly differ from those of the existing theoretical models and similarity solutions for axisymmetric currents such as the one derived for a fixed point source in Appendix D. In addition, these solutions are only valid when the plume head aspect ratio |$\mathcal {H}/\mathcal {R}$| is much less than 1, that is far downstream the source. In our simulations, upon its arrival beneath the plate, the plume head aspect ratios |$\mathcal {H}/\mathcal {R}$| are |$\mathcal {O}(1)$| as detailed in Table 2 and by the time the plume reaches the slab, |$\mathcal {H}/\mathcal {R}$| has decreased to 0.1–0.3 only. So, regardless of the external mantle flow, we cannot expect that plume heads follow a propagation as t3/5 in the early times. The existence of a subsequent time window (before slab/plume impact) that is characterized a propagation as t3/5 depends on whether the initial slab–plume distance is sufficiently long and the slab retreat rate is sufficiently small compared to that of plume spreading. The importance of the plume head aspect ratio |$\mathcal {H}/\mathcal {R}$| is shown in Fig. 17 by the plume-only runs P2 and P5, which are free from any background mantle flow and thus exhibit axisymmetric plumes. Plume-only run P2 does not converge to a propagation as t3/5 until about the last 2/5th of the simulated propagation. For plume-only run P5, it is even less, until about the last 1/5th of the simulated propagation. The comparison between slab–plume runs 2 and 5 and their plume-only equivalent runs P2 and P5 further shows the increase in lateral spreading in the slab–plume simulations comes at the expense of the width growth. In run 5 featuring a ‘weak plume’, the head of the plume is affected as early as during its rise. When the head reaches the plate (t ≈ 650 s), the plume half-width |$w(t)/2$| differs by 8 % from its analogue without slab (Run P5). The slab’s influence thus took place while the slab was between 18 and 11 cm away from the plume source xp (i.e. when 23 < xh < 30 cm as shown in Fig. 5). We note that the slab’s influence also affected the conduit during the plume ascent by tilting it as shown in Figs 14 and A2. Figure 16. Open in new tabDownload slide Asymmetry of the propagated plume head relative to the source position before its subduction. Top views for simulations representative of the range of buoyancy ratios Bs/Bp = [3, 30] (Bs/Bp decreases from the top to the bottom panels; see Table 2). The green circles indicate the location of the conduit beneath. The green line locates the maximum extent of the plume in the z direction. The insets show a side view of the conduit from the back wall. Figure 16. Open in new tabDownload slide Asymmetry of the propagated plume head relative to the source position before its subduction. Top views for simulations representative of the range of buoyancy ratios Bs/Bp = [3, 30] (Bs/Bp decreases from the top to the bottom panels; see Table 2). The green circles indicate the location of the conduit beneath. The green line locates the maximum extent of the plume in the z direction. The insets show a side view of the conduit from the back wall. Figure 17. Open in new tabDownload slide (a) Half-width |$w(t)/2=(x_u(t)-x_d(t))/2$| as a function of time t in slab/plume and plume-only runs 2, 5 and P2 and P5, respectively. (b) |$z_l(t)$| as a function of time in slab/plume and plume-only runs 2, 5 and P2 and P5, respectively. The dotted black lines are the slopes for an axisymmetric spreading. Figure 17. Open in new tabDownload slide (a) Half-width |$w(t)/2=(x_u(t)-x_d(t))/2$| as a function of time t in slab/plume and plume-only runs 2, 5 and P2 and P5, respectively. (b) |$z_l(t)$| as a function of time in slab/plume and plume-only runs 2, 5 and P2 and P5, respectively. The dotted black lines are the slopes for an axisymmetric spreading. 5 DISCUSSION In this study, we examined the dynamics of head-to-toe mantle plumes at the rear of a retreating slab using 3-D numerical simulations. Nine slab–plume simulations were performed over a range of slab–plume buoyancy flux ratios. The influence of the plume on the slab was found to be a long-distance influence causing the decrease of the retreat rate when plume and slab are as far apart as a few hundreds of kilometres (400–850 km) as identified by the start of the slab’s retreat phase III in Section 3 T scaled to nature (× 50 km). Similarly, the slab’s influence on the plume was found to be long-distance as exemplified in run 5 where the slab was at least 550 km apart from the plume source when it influenced the plume dynamics as shown in Section 4.3. We suspect that, in our simulations, these latter length scales of influence actually underestimate the maximum plume–slab distance of interaction due to the limitation imposed by the plate length and consequently the initial plume source-slab hinge distance. We suspect this as run 7 in the long-plate configuration had a length of influence 36 % longer than that of run 2 in the short-plate configuration with a similar buoyancy as pointed out in Section 3. These results are interesting because one could have assumed that the maximum slab-plume interaction distance is primarily dictated by the mantle depth below the plate ( H − d =335 km) controlling the convection cell, in which case the plume influence on the slab would not have started before plume and slab are 335 km apart. Finding a scaling for the maximum plume–slab interaction distance remains a question to be clarified in the future. These simulations compare well with previous laboratory experiments (Mériaux et al. 2015). In particular, the two modelling approaches show a similar impact of the head once it has reached the slab. This causes an additional decrease of the retreat rate over a range of slab–plume buoyancy flux ratios (Bs/Bp⪅ 25). However, by the virtue of having a complete space–time record of the velocity field, numerical simulations made it possible to better interpret the dynamics of the slab–plume interaction and more specifically the influence of the feeding conduit on the spreading plume head and vice versa. We thus showed that the plume conduit trajectory at the rear of a subduction in retreat mode is a result of both the slab– and plume(head)–induced mantle flow. The plume–induced flow could contribute to a 30 % increase in the shear affecting the conduit in our simulations. To our knowledge, this plume auto-shearing mechanism increasing conduit deflection has never been reported. Plume conduits in our simulations happened to be strongly tilted to up to 55° relative to the vertical, which is known to favour the development of Rayleigh-Taylor instability. For instance, according to Lister et al. (2011), the wavelength of instability along a buoyant conduit inclined at a 40° angle rising through a quiet and unbounded ambient would be ≈4 cm at our ambient/plume viscosity contrast. This wavelength is less than the 10 cm conduit length over the depth of mantle below the plate (6.7 cm) at that inclination. Yet, we do not observe any instabilities. It might be that the ambient flow suppresses them, but this needs to be further investigated. We also demonstrate that the conduit trajectory can be approximated by a third-order polynomial function of depth that is derived from the quadratic slab-only-induced poloidal flow profile when Bs/Bp ≤ 9.5. In such cases, we provide a rule of thumb to estimate the offset of the conduit nearby the base of the plume, which is given by eq. (6). Ultimately, the plume spreads while being fed by a conduit moving away from the slab. This motion increases as the slab–plume buoyancy flux ratio increases. In our simulations, the conduit beneath the head moved by ≈ 2–5 cm in ∼1000 s, that is ≈ 100–250 km in ∼30 Myr. In addition to the tilt, the plume conduit gradually flares as a function of height. In the case of strong plumes, the deformation can be associated with exponential stretching due to the toroidal flow. The conduit tilt and flare increase during slab retreat, changing the feeding conditions to the plume head. In particular, they cause an increasingly unbalanced azimuthal distribution of the plume source flux into the head. The plume flux released into the head by the conduit progressively flows in the direction parallel to the trench and away from the slab. The combined effect of the feeding conditions and ambient flow in the mantle results in an asymmetric spreading of the head essentially in the direction parallel to the trench. Based on depth velocity profiles, the effect of the ambient flow on the head is however expected to be less significant than its effect on the conduit in our simulations. The offset of the conduit with respect to the plume head and the propagation of the plume head parallel to the trench are features to be expected in nature. In the latter context, the event envisioned by geodynamics studies between the Yellowstone plume and the retreating Farallon slab is quite unique as the rising plume is thought to have directly impacted the apex of the slab without prior propagation beneath the lithosphere (e.g. Leonard & Liu 2016). Less conjecturable is a process of plume-modified orogenesis that begins by a stage where the plume impinges the oceanic lithosphere and propagates towards the retreating slab before tearing it from beneath once subducted, as proposed in the geologic past for the Abitibi–Wawa province (Wyman et al. 2002), or the southern Central Andes (Gianni et al. 2017). Today, there is no clear example of such a stage of slab–plume interaction. For instance, the plume heads of the San Felix and Juan Fernandez hotspots in South America have already been subducted beneath the Nazca slab. However, our simulations show that the dominant propagation of the plume head in the direction parallel to the trench prior to its impact with the slab implies that a front of buoyant material could disturb subduction. In our simulations, the trench-parallel plume front represents the widening of the plume head by two to two and a half times since the start of spreading. The front extent |$(\sim 2z_l)$| ranged from 10 to 20 cm, that is 500–1000 km, at respective slab-plume buoyancy ratios 30 and 3. These bounds for the extent of the trench-parallel plume front are likely underestimates compared to those that would result from the interaction in Nature because, even if they will be changes to the rate of propagation of the plume head due to the presence of the asthenosphere and shearing of the plume by the plate velocity, the initial thermal plume heads will be larger than their compositional counterparts (∼800–1200 km) due to thermal entrainment during the ascent (Griffiths & Campbell 1990). The slab–plume interaction of a starting plume at the rear of a retreating slab is thus bound to be very different from the slab–plume interaction with the only subsequent products of the plume conduit (e.g. a narrower ridge). The initial impact of the plume head with the slab and the plume’s subduction may even more so contribute to slab tearing or breaking, processes that are often invoked in the subsequent stages of plume-modified orogenies. In the context of this study, we also showed that the application of gravity current similarity solutions for the spreading of the plume head is only valid once the plume head characteristic length |$\mathcal {R}$| is at least one order of magnitude less that the characteristic head thickness |$\mathcal {H}$|⁠. In our simulations, the convergence to the similarity solutions required that R ≈ 2R0, where R0 is the initial plume radius impinging the lithosphere, which took ∼800 s, that is ∼23 Myr. So, a typical plume head of radius 500 km impinging the lithosphere would need to propagate until its head has reached a radius of 1000 km before the similarity solutions could be applied. It is not clear if this transient would be shorter when the plume head spreads in a mantle layer of viscosity similar to that of the plume (μp ∼ μm ; e.g. asthenosphere) rather than a mantle layer of larger viscosity (μp ≪ μm; e.g. upper mantle) and with an imposed velocity at the top surface of the layer. However, in the absence of any background motion, when the similarity solutions for an axisymmetric propagation apply (i.e. |$\mathcal {H}/\mathcal {R}\ll 1$|⁠), a slower decrease of the aspect ratio is predicted when μp ∼ μm. The aspect ratio |$\mathcal {H}/\mathcal {R}$| decreases as t−1/2 when μp ∼ μm, whereas it varies as t−4/5 when μp ≪ μm (see Appendix D). In any case, the early stages of the plume head spreading are controlled by the near-source flow conditions given by the instantaneous flux from the source. We hope that this study can be an incentive for further investigations on the conduit supply conditions to a plume head or plume conduit material during interaction with the lithosphere as they are likely to impact the topography of hotspot swells and melting. ACKNOWLEDGEMENTS The numerical work was supported by computational resources provided in part by the Australian Government through the use of the HPC Raijin cluster of the National Computational Infrastructure (NCI) under the National Computational Merit Allocation Scheme (NCMAS) project qk0, and in part by the Monash eResearch Centre and eSolutions-Research Support Services through the use of the MonARCH HPC Cluster. Author D. May acknowledges financial support from the European Research Council under the European Union’s Seventh Framework Programme(FP7/2007-2013)/ERC Grant Agreement Number 279925 and the Alfred P. Sloan Foundation through the Deep Carbon Observatory (DCO) Modeling and Visualization Forum. We thank Editor Professor Juan Carlos Afonso and three referees for their feedbacks, which contributed to the improvement of the manuscript. REFERENCES Balay S. , Gropp W.D., McInnes L.C., Smith B.F., 1997 . Efficient management of parallelism in object oriented numerical software libraries , in Modern Software Tools in Scientific Computing , pp. 163 – 202 ., Birkhäuser Press . 10.1007/978-1-4612-1986-6_8 Google Scholar Crossref Search ADS Google Scholar Google Preview WorldCat COPAC Crossref Balay S. et al. ., 2017a . PETSc users manual, Tech. Rep. ANL-95/11 - Revision 3.8 , Argonne National Laboratory . Balay S. et al. ., 2017b . PETSc Web page , http://www.mcs.anl.gov/petsc . OpenURL Placeholder Text WorldCat Betts P.G. et al. ., 2009 . Mesoproterozoic plume-modified orogenesis in eastern Precambrian Australia , Tectonics , 28 ( 3 ), doi:10.1029/2008TC002325.doi: 10.1029/2008TC002325 OpenURL Placeholder Text WorldCat Crossref Betts P.G. , Mason W.G., Moresi L., 2012 . The influence of a mantle plume head on the dynamics of a retreating subduction zone , Geology , 40 ( 8 ), 739 – 742 . 10.1130/G32909.1 Google Scholar Crossref Search ADS WorldCat Crossref Brendan Murphy J. , Oppliger G.L., Brimhall G.H.Jr, Hynes A., 1998 . Plume-modified orogeny: an example from the western United States , Geology , 26 ( 8 ), 731 – 734 . 10.1130/0091-7613(1998)026&lt;0731:PMOAEF&gt;2.3.CO;2 Google Scholar Crossref Search ADS WorldCat Crossref Chang S.-J. , Ferreira A.M., Faccenda M., 2016 . Upper-and mid-mantle interaction between the Samoan plume and the Tonga–Kermadec slabs , Nat. Commun. , 7 , 10799 . 10.1038/ncomms10799 Google Scholar Crossref Search ADS PubMed WorldCat Crossref Crameri F. et al. ., 2012 . A comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the “sticky air” method , Geophys. J. Int. , 189 ( 1 ), 38 – 54 . 10.1111/j.1365-246X.2012.05388.x Google Scholar Crossref Search ADS WorldCat Crossref Feighner M.A. , Richards M.A., 1995 . The fluid dynamics of plume-ridge and plume-plate interactions: an experimental investigation , Earth planet. Sci. Lett. , 129 ( 1–4 ), 171 – 182 . 10.1016/0012-821X(94)00247-V Google Scholar Crossref Search ADS WorldCat Crossref Gianni G.M. , García H.P., Lupari M., Pesce A., Folguera A., 2017 . Plume overriding triggers shallow subduction and orogeny in the southern central Andes , Gondw. Res. , 49 , 387 – 395 . 10.1016/j.gr.2017.06.011 Google Scholar Crossref Search ADS WorldCat Crossref Griffiths R.W. , Campbell I.H., 1990 . Stirring and structure in mantle starting plumes , Earth planet. Sci. Lett. , 99 ( 1–2 ), 66 – 78 . 10.1016/0012-821X(90)90071-5 Google Scholar Crossref Search ADS WorldCat Crossref Huppert H.E. , 1982 . The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface , J. Fluid Mech. , 121 , 43 – 58 . 10.1017/S0022112082001797 Google Scholar Crossref Search ADS WorldCat Crossref Kaluza O. , Moresi L., Mansour J., Barnes D.G., 2019 . Okaluza/lavavu: v1.3.2 . 10.1093/gji/ggaa222 Crossref Kerr R.C. , Mériaux C., 2004 . Structure and dynamics of sheared mantle plumes , Geochem. Geophys. Geosyst. , 5 ( 12 ), doi:10.1029/2004GC000749. 10.1029/2004GC000749 OpenURL Placeholder Text WorldCat Crossref Lee C. , Lim C., 2014 . Short-term and localized slab-plume interaction explains the genesis of Abukuma Adakite in northeastern Japan , Earth planet. Sci. Lett. , 396 , 116 – 124 . 10.1016/j.epsl.2014.04.009 Google Scholar Crossref Search ADS WorldCat Crossref Leonard T. , Liu L., 2016 . The role of a mantle plume in the formation of yellowstone volcanism , Geophys. Res. Lett. , 43 ( 3 ), 1132 – 1139 . 10.1002/2015GL067131 Google Scholar Crossref Search ADS WorldCat Crossref Lister J.R. , Kerr R.C., 1989 . The propagation of two-dimensional and axisymmetric viscous gravity currents at a fluid interface , J. Fluid Mech. , 203 , 215 – 249 . DOI: 10.1017/S0022112089001448 Google Scholar Crossref Search ADS WorldCat Crossref Lister J.R. , Kerr R.C., Russell N.J., Crosby A., 2011 . Rayleigh–Taylor instability of an inclined buoyant viscous cylinder , J. Fluid Mech. , 671 , 313 – 338 . 10.1017/S0022112010005689 Google Scholar Crossref Search ADS WorldCat Crossref Manga M. , 1996 . Mixing of heterogeneities in the mantle: effect of viscosity differences , Geophys. Res. Lett. , 23 ( 4 ), 403 – 406 . 10.1029/96GL00242 Google Scholar Crossref Search ADS WorldCat Crossref Mansour J. et al. ., 2020 . Underworld2: Python geodynamics modelling for desktop, hpc and cloud , J. Open Sour. Soft. , 5 ( 47 ), 1797 . 10.21105/joss.01797 Google Scholar Crossref Search ADS WorldCat Crossref Mériaux C. , Jaupart C., 1995 . Simple fluid dynamic models of volcanic rift zones , Earth planet. Sci. Lett. , 136 ( 3–4 ), 223 – 240 . Google Scholar Crossref Search ADS WorldCat Mériaux C. , Mériaux A.-S., Schellart W., Duarte J., Duarte S., Chen Z., 2016 . Mantle plumes in the vicinity of subduction zones , Earth planet. Sci. Lett. , 454 , 166 – 177 .. 10.1016/j.epsl.2016.09.001 Google Scholar Crossref Search ADS WorldCat Crossref Mériaux C.A. , Duarte J.C., Schellart W.P., Mériaux A.-S., 2015 . A two-way interaction between the Hainan plume and the manila subduction zone , Geophys. Res. Lett. , 42 ( 14 ), 5796 – 5802 . 10.1002/2015GL064313 Google Scholar Crossref Search ADS WorldCat Crossref Mériaux C.A. , May D.A., Mansour J., Chen Z., Kaluza O., 2018 . Benchmark of three-dimensional numerical models of subduction against a laboratory experiment , Phys. Earth planet. Inter. , 283 , 110 – 121 . 10.1016/j.pepi.2018.07.009 Google Scholar Crossref Search ADS WorldCat Crossref Moresi L. , Quenette S., Lemiale V., Mériaux C., Appelbe B., Mühlhaus H.-B., 2007 . Computational approaches to studying non-linear dynamics of the crust and mantle , Phys. Earth planet. Inter. , 163 ( 1 ), 69 – 82 . 10.1016/j.pepi.2007.06.009 Google Scholar Crossref Search ADS WorldCat Crossref Moresi L. , Mansour J., Giordani J., Farrington R., Beucher R., 2018 . Underworld 2, Note: future releases will automatically be uploaded via github. This is a manual upload of the 2018 Q1 release. The master URL for underworld is 10.5281/zenodo.1436039 . 10.1093/gji/ggaa222 Crossref Olson P. , 1990 . Hot spots, swells and mantle plumes , in Magma Transport and Storage , Ryan M.P., pp. 33 – 51 ., Wiley . OpenURL Placeholder Text WorldCat Ribe N. , Christensen U., 1994 . Three-dimensional modeling of plume-lithosphere interaction , J. geophys. Res. , 99 ( B1 ), 669 – 682 . 10.1029/93JB02386 Google Scholar Crossref Search ADS WorldCat Crossref Richards M.A. , Griffiths R.W., 1988 . Deflection of plumes by mantle shear flow: experimental results and a simple theory , Geophys. J. Int. , 94 ( 3 ), 367 – 376 . 10.1111/j.1365-246X.1988.tb02260.x Google Scholar Crossref Search ADS WorldCat Crossref Schmeling H. et al. ., 2008 . A benchmark comparison of spontaneous subduction models?towards a free surface , Phys. Earth planet. Inter. , 171 ( 1 ), 198 – 223 . 10.1016/j.pepi.2008.06.028 Google Scholar Crossref Search ADS WorldCat Crossref Sleep N.H. , 1996 . Lateral flow of hot plume material ponded at sublithospheric depths , J. geophys. Res. , 101 ( B12 ), 28065 – 28083 . 10.1029/96JB02463 Google Scholar Crossref Search ADS WorldCat Crossref Spence D. , Ockendon J., Wilmott P., Turcotte D., Kellogg L., 1988 . Convective mixing in the mantle: the role of viscosity differences , Geophys. J. Int. , 95 ( 1 ), 79 – 86 . 10.1111/j.1365-246X.1988.tb00452.x Google Scholar Crossref Search ADS WorldCat Crossref Tassara S. , et al. ., 2017 . Plume-subduction interaction forms large auriferous provinces , Nat. Commun. , 8 ( 1 ), 843 . 10.1038/s41467-017-00821-z Google Scholar Crossref Search ADS PubMed WorldCat Crossref Taylor G.I. , 1934 . The formation of emulsions in definable fields of flow , Proc. R. Soc. Lond., A , 146 ( 858 ), 501 – 523 . . Google Scholar Crossref Search ADS WorldCat Wyman D. , Kerrich R., Polat A., 2002 . Assembly of archean cratonic mantle lithosphere and crust: plume–arc interaction in the Abitibi–Wawa subduction–accretion complex , Precambrian Res. , 115 ( 1-4 ), 37 – 62 . 10.1016/S0301-9268(02)00005-0 Google Scholar Crossref Search ADS WorldCat Crossref APPENDIX A: TIME SEQUENCE OF PLUME/SLAB RUNS Here, we present two sets of snapshots of the system at times when the plume arrives beneath the plate yp = 6 cm and then when plume(head)-slab distances are xh − xu = 8, 4, 1, −0.5 cm. The same plume(head)-slab distances are shown in Figs A1 and A2 so that they can be compared. Note that the evolution is similar to that shown in previous laboratory studies (e.g. Mériaux et al. 2015, 2016). The two sets were generated from runs 2 and 5, which represent two-end members of a strong and weak plume relative to the slab. Figure A1. Open in new tabDownload slide Slab–plume evolution in run 2 seen from two side views and a top view. Figure A1. Open in new tabDownload slide Slab–plume evolution in run 2 seen from two side views and a top view. Figure A2. Open in new tabDownload slide Slab–plume evolution in run 5 seen from two side views and a top view. Figure A2. Open in new tabDownload slide Slab–plume evolution in run 5 seen from two side views and a top view. APPENDIX B: STEADY PLUME CONDUIT TRAJECTORY IN A PARABOLIC SHEAR PROFILE The trajectory of a plume conduit results from the combination of the plume rise velocity and the background mantle flow. For instance, plume conduits rising from a fixed source and bent over by a constant background mantle flow following a depth-dependent linear velocity profile are expected to reach a steady parabolic trajectory (Richards & Griffiths 1988). As shown in Fig. 7, the background mantle flow in the x (negative) direction at the rear of the slab has a parabolic profile in the mantle. The velocity vx(y) can be approximated by a polynomial function of order 2 over nearly the entire depth of mantle fluid Hm given by $$\begin{eqnarray*} v_x(y)= v_0 \left(\frac{y^2}{H_m^2} - b \frac{y}{H_m} \right), \end{eqnarray*}$$(B1) where |$v_0$| is a velocity pre-factor and b is a positive dimensionless constant. The velocity pre-factor |$v_0$| is related to the maximum shearing velocity in the layer |$v_{max}$| by |$v_0=-4v_{max}/b$|⁠. In our simulations, b was found to be close to 1. Assuming that each plume conduit element rises at a velocity that is the vector sum of the background velocity vx(y) and the Stokes velocity |$v_S=k \Delta \rho g r_p^2/ \mu _m$|⁠, where k = 0.54 (Richards & Griffiths 1988), the conduit trajectory can be defined as $$\begin{eqnarray*} x(y) = \int _0^y \frac{v_x(y)}{v_S} dy. \end{eqnarray*}$$(B2) Upon integration of eq. (B2) with eq. (B1), x(y) takes the form $$\begin{eqnarray*} x(y) = \frac{v_0}{v_S} \left(\frac{y^3}{3H_m^2} - b \frac{y^2}{2H_m} \right). \end{eqnarray*}$$(B3) Fig. B1 shows the conduit trajectory generated in the slab-induced ambient shear predicted by the application of eqs (B1) and (B3). Note that the final offset of the conduit is $$\begin{eqnarray*} x(H_m)=\frac{v_0H_m}{3v_S} \left(1 - \frac{3b}{2} \right). \end{eqnarray*}$$(B4) Figure B1. Open in new tabDownload slide (a) Profile of the horizontal velocity |$v_x$| at a fixed position z = 0 cm and x = xp + 6 = 18 cm at the rear of the slab as function of height y for t ∼700 s in slab–only run 1 (solid black line). Eq. (B1), where Hm = H − d and |$v_0$| and b have been fitted to the velocity profile is the dashed magenta line. (b) Associated shape of the conduit calculated with eq. (B3) for four plume radii. Figure B1. Open in new tabDownload slide (a) Profile of the horizontal velocity |$v_x$| at a fixed position z = 0 cm and x = xp + 6 = 18 cm at the rear of the slab as function of height y for t ∼700 s in slab–only run 1 (solid black line). Eq. (B1), where Hm = H − d and |$v_0$| and b have been fitted to the velocity profile is the dashed magenta line. (b) Associated shape of the conduit calculated with eq. (B3) for four plume radii. APPENDIX C: CONDUIT EVOLUTION FROM AXISYMMETRIC TO ELLIPTIC In the horizontal x−z planes, the deformation of the conduit resembles that of a viscous circular blob into an ellipsoid while in pure shear (i.e. |$v_x=\dot{\varepsilon } x$| and |$v_z=-\dot{\varepsilon } z$| as (x2 + z2)1/2 → ∞, where |$\dot{\varepsilon }$| is the strain rate and |$v_x$| and |$v_z$| are the velocity components in the x and z directions). The difference is that this incoming flow is not symmetric as in pure shear. As shown in Fig. C1, the deformation in pure shear has first been described by Taylor (1934) and then by Spence et al. (1988) and Manga (1996) in the context of 2-D mixing in the convective mantle. These studies showed that a pure shear flow leads to exponential stretching of a viscous blob following: $$\begin{eqnarray*} \mbox{ln}\left( \frac{r_c^z(t)}{r_c^z(0)} \right) + \frac{ \lambda r_c^x(0)}{r_c^z(0)}\left( 1-\frac{(r_c^z(0))^2}{ (r_c^z(t))^2}\right) =\dot{\varepsilon } t, \end{eqnarray*}$$(C1) where λ is the viscosity ratio between the blob and the surrounding mantle, and |$r_c^z$| and |$r_c^x$| are the semi-axes of the ellipse. According to eq. (C1), the amount of stretching depends on the initial aspect ratio |$r_c^x(0)/r_c^z(0)$| and the viscosity contrast λ. In the limit where |$\lambda r_c^x(0)/ r_c^z(0) \ll 1$|⁠, eq. (C1) becomes $$\begin{eqnarray*} r_c^z(t)= r_c^z(0) \mbox{exp}\left(\dot{\varepsilon } t \right). \end{eqnarray*}$$(C2) Figure C1. Open in new tabDownload slide 2-D stretching in pure shear flow. Figure C1. Open in new tabDownload slide 2-D stretching in pure shear flow. In this study, λ = μp/μm = 0.0455, so eq. (C2) is valid when |$r_c^x(0)/ r_c^z(0) \ll 22$|⁠, which is satisfied for the conduit and the plume head at all times in our study. Assuming that the initial circular conduit is deformed into an ellipse, the tracking of length |$r_c^z$| over time at different head heights is sufficient for inferring the evolution of the conduit |$r_c^x$| in the x direction since the volume flux must be conserved independently of height. Whereas the laminar flux through an axisymmetric conduit follows eq. (4), the same flux Q through an elliptical conduit is given by $$\begin{eqnarray*} Q=\frac{ \pi \Delta \rho _p g}{4\mu _p} \frac{(r_c^x)^3 + (r_c^z)^3}{(r_c^x)^2 + (r_c^x)^2}, \end{eqnarray*}$$(C3) where |$r_c^x$| and |$r_c^z$| are the semi axes in the x and z directions, respectively (Mériaux & Jaupart 1995). Applying the conservation of the flux thus thus gives the relationship $$\begin{eqnarray*} \frac{r_p^4}{2}= \frac{(r_c^x)^3 (r_c^z)^3}{(r_c^x)^2 + (r_c^z)^2}. \end{eqnarray*}$$(C4) Eq. (C4) is a third-order polynomial function of |$r_c^x$| following: $$\begin{eqnarray*} (r_c^x)^3 - \frac{r_p^4}{2(r_c^z)^3} (r_c^x)^2 - \frac{r_p^4}{2r_c^z}= 0. \end{eqnarray*}$$(C5) The discriminant Δ of eq. (C5) given by $$\begin{eqnarray*} \Delta = \frac{17}{243} \frac{r_p^{24}}{(r_c^z)^{18}} + \frac{4r_p^{16}}{27 (r_c^z)^{10} } \end{eqnarray*}$$(C6) is positive. Therefore, there is only one real root of eq. (C5) $$\begin{eqnarray*} r_c^x= \left(-\frac{X}{2} + \sqrt{\Delta } \right) ^{1/3} + \left(-\frac{X}{2} - \sqrt{\Delta } \right) ^{1/3} + \frac{r_p^4}{6(r_c^z)^3}, \end{eqnarray*}$$(C7) where $$\begin{eqnarray*} X= \frac{2r_p^4}{54(r_c^z)^3} - \frac{r_p^4}{2r_c^z}. \end{eqnarray*}$$(C8) APPENDIX D: PROPAGATION OF VISCOUS GRAVITY CURRENTS Let’s first consider a system where a fixed point source, delivering at volumetric rate Qtα, a fluid of density ρ − Δρ and viscosity μi, between a fluid of density ρ and viscosity μa ≫ μi, and a rigid horizontal surface (see Fig. D1a). In the absence of any background motion in the ambient fluid, the intruding fluid spreads axisymmetrically and the viscous dissipation primarily occurs in the ambient fluid. Assuming that the scale for thickness of the current is much less than the length scale of the current (⁠|$\mathcal {H}\ll \mathcal {R}$|⁠), a scaling law for the propagation of the gravity current can be derived from a force balance between the gravity and viscous forces and the volume conservation. Figure D1. Open in new tabDownload slide Axisymmetric viscous gravity current fed by a fixed point source. View in any vertical planes. Figure D1. Open in new tabDownload slide Axisymmetric viscous gravity current fed by a fixed point source. View in any vertical planes. The horizontal gravity and viscous forces Fg and Fv are, respectively, given by $$\begin{eqnarray*} F_g \sim \Delta \rho g h^2 r_n, \end{eqnarray*}$$(D1) where h(t) is the current thickness and rn(t) is the radius of the current and $$\begin{eqnarray*} F_v \sim \mu _a \frac{r_n}{t} r_n. \end{eqnarray*}$$(D2) Balancing eqs (D1) and (D2) results in $$\begin{eqnarray*} r_n \sim \frac{\Delta \rho g h^2}{\mu _a} t. \end{eqnarray*}$$(D3) When eq. (D3) is combined with a global volume constraint $$\begin{eqnarray*} r_n^2 h \sim Qt^{\alpha }, \end{eqnarray*}$$(D4) a solution for the radius of the current is found. When α = 0, a fixed volume is released. When α = 1, the current is fed at a constant rate. It follows: $$\begin{eqnarray*} r_n(t) \sim \left(\frac{\Delta \rho g Q^2}{\mu _a} \right)^{1/5} t^{(2\alpha +1)/5}, \end{eqnarray*}$$(D5) and $$\begin{eqnarray*} h(t) \sim \left(\frac{\Delta \rho g }{\mu _a Q^{1/2}} \right)^{-2/5} t^{(\alpha -2)/5}. \end{eqnarray*}$$(D6) This solution can also be derived using a lubrication-theory approximation and similarity solutions as shown by Lister & Kerr (1989). In laboratory experiments of axisymmetric currents at fixed volume flux, Feighner & Richards (1995) found that $$\begin{eqnarray*} r_n(t) = 0.79 \left(\frac{\Delta \rho g Q^2}{\mu _a} \right)^{1/5} t^{3/5}. \end{eqnarray*}$$(D7) Note that in this case, h/rn also evolves as ∝t−4/5. Eventually, as the current itself lengthens, the viscous dissipation inside the current increases and another regime of propagation is expected [see for instance (Huppert 1982; Kerr & Mériaux 2004)]. In the case, where the dissipation occurs within an axisymmetric current, h/rn evolves as ∝t−1/2. © The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Mantle plume dynamics at the rear of a retreating slab
JF - Geophysical Journal International
DO - 10.1093/gji/ggaa222
DA - 2020-08-01
UR - https://www.deepdyve.com/lp/oxford-university-press/mantle-plume-dynamics-at-the-rear-of-a-retreating-slab-blc6Yg2agM
SP - 1146
VL - 222
IS - 2
DP - DeepDyve
ER -