The mantle wedge's transient 3‐D flow regime and thermal structure

Arc volcanism, volatile cycling, mineralization, and continental crust formation are likely regulated by the mantle wedge's flow regime and thermal structure. Wedge flow is often assumed to follow a regular corner‐flow pattern. However, studies that incorporate a hydrated rheology and thermal buoyancy predict internal small‐scale‐convection (SSC). Here, we systematically explore mantle‐wedge dynamics in 3‐D simulations. We find that longitudinal “Richter‐rolls” of SSC (with trench‐perpendicular axes) commonly occur if wedge hydration reduces viscosities to ≲1·1019 Pa s, although transient transverse rolls (with trench‐parallel axes) can dominate at viscosities of ∼5·1018−1·1019 Pa s. Rolls below the arc and back arc differ. Subarc rolls have similar trench‐parallel and trench‐perpendicular dimensions of 100–150 km and evolve on a 1–5 Myr time‐scale. Subback‐arc instabilities, on the other hand, coalesce into elongated sheets, usually with a preferential trench‐perpendicular alignment, display a wavelength of 150–400 km and vary on a 5–10 Myr time scale. The modulating influence of subback‐arc ridges on the subarc system increases with stronger wedge hydration, higher subduction velocity, and thicker upper plates. We find that trench‐parallel averages of wedge velocities and temperature are consistent with those predicted in 2‐D models. However, lithospheric thinning through SSC is somewhat enhanced in 3‐D, thus expanding hydrous melting regions and shifting dehydration boundaries. Subarc Richter‐rolls generate time‐dependent trench‐parallel temperature variations of up to ∼150 K, which exceed the transient 50–100 K variations predicted in 2‐D and may contribute to arc‐volcano spacing and the variable seismic velocity structures imaged beneath some arcs.


Introduction
The majority of Earth's seismicity and volcanism occurs at destructive plate margins, where subducting plates (slabs) sink below overriding (upper) plates into Earth's mantle, transporting volatiles and other elements to depth [e.g., Stern, 2002]. In the underlying mantle wedge, viscous drag, induced by the subducting slab, forces mantle material to flow from the back-arc region toward the wedge corner, where downgoing and upper plates meet [e.g., Davies and Stevenson, 1992]. The resultant high wedge-corner temperatures lead to mineral dehydration in the downgoing slab, delivering water to the overlying mantle, thus facilitating melt formation and magmatism [e.g., Gill, 1981;Davies and Stevenson, 1992;Tatsumi and Eggins, 1995;Stern, 2002;van Keken et al., 2002;Wilson et al., 2014]. Numerical models that simulate the mantle wedge's flow regime and thermal structure aim to reproduce these conditions [e.g., Davies and Stevenson, 1992;van Keken et al., 2002;Syracuse et al., 2010;Le Voci et al., 2014;Wilson et al., 2014], with the goal of better understanding the key controls on the location and style of arc volcanism and the observed variability between different subduction zones.
Previous numerical studies have generally focussed on 2-D simulations of mantle wedge flow, kinematically driven by the subducting plate, with most neglecting the role of a hydrated wedge rheology and local thermal buoyancy. Such models yield wedge and slab surface temperatures that agree with many observations [e.g., Kneller et al., 2007;Plank et al., 2009;Syracuse et al., 2010;Long and Becker, 2010;Hebert and Gurnis, 2010]. However, they cannot explain the complex wedge seismic structure imaged beneath some volcanic arcs, temporal and trench-parallel variability in arc volcanism, or the thin lithosphere and high heat-flow observed in many back-arc regions [e.g., Tamura et al., 2002;Schurr et al., 2003;Currie and Hyndman, 2006].

Key Points:
Small-scale-convection (SSC) occurs in the mantle wedge at viscosities below $1.10 19 Pa s SSC thins the lithosphere, expands hydrous melting regions, and shifts dehydration boundaries The resulting 150 K trench-parallel temperature variations may contribute to arc-volcano spacing Supporting Information: Figure S1  Under hydrous conditions, viscosities are likely substantially lower than those of dry mantle [e.g., Karato and Wu, 1993;Hirth and Kohlstedt, 1996], leading to small-scale convection (SSC) that is driven by gravitational instabilities from the base of the upper plate [e.g., Honda and Saito, 2003;Wirth and Korenaga, 2012;Le Voci et al., 2014]. Although some recent studies have challenged the rheological influence of water in the wedge [e.g., Fei et al., 2013;Girard et al., 2013], and furthermore, melting extracts water and, hence, may dry the wedge [e.g., Hebert et al., 2009], low seismic velocity anomalies, low seismic attenuation, and the back-arc surface topography at a number of subduction zones are most easily explained with a low-viscosity mantle wedge, potentially extending hundreds of kilometers from the trench below the back arc [e.g., Billen and Gurnis, 2001;Currie and Hyndman, 2006;Wiens et al., 2008;Greve et al., 2014].
In a previous 2-D study [Le Voci et al., 2014], we systematically examined the effect of viscosity (as a function of wedge hydration), subduction velocity, slab dip, and upper plate age (thickness) on the mantle wedge's flow regime and thermal structure. Consistent with previous work [e.g., Honda and Yoshida, 2005;Honda et al., 2010;Wirth and Korenaga, 2012], we found SSC to be a common occurrence when wedge hydration lowered mantle viscosities to below Շ5 Á 10 18 Pa s. In addition, such low viscosities needed to extend over lateral distances that exceeded instability wavelengths (i.e., over distances greater than $150 km). When present, SSC led to substantial lithospheric thinning and a transient 50-100 K variability in wedge temperatures, which is sufficient to affect melting and dehydration. While net back-arc lithospheric thinning was largest when slab velocities were highest, drips were most pronounced when subducting plate velocities and, thereby, the shearing of drips by background corner-flow, was lowest.
As in other kinematically driven wedge models, thinning of the upper plate lithosphere was strongest in the wedge corner, more or less above where the slab reaches a depth that in observations corresponds to volcanic-arc locations [e.g., England et al., 2004]. Although models like ours do not predict arc locations, we will refer to the region as the arc/subarc region. Dislocation-creep enhanced strain-rates and elevated temperatures in the wedge corner promotes the development of a ''pinch-zone,'' where isotherms are compressed against the slab's surface [e.g., Kincaid and Sacks, 1997]. This is also a zone of differential thinning between the subarc and subback-arc regions, as the mantle wedge penetrates upwards as a narrow hot tongue into the overriding plate, resulting in a local decrease in upper plate thickness. The size of this zone is affected by the implementation of the decoupled slab/upper plate interface [e.g., Conder, 2005;Arcay et al., 2008], which previous studies indicate must extend to $80 km depth in order to satisfy observations of low surface heat flow and low seismic attenuation in the fore-arc region [e.g., Wada and Wang, 2009;Syracuse et al., 2010]. In our 2-D models, the pinch zone was most pronounced for cases with the highest wedge viscosities, subduction velocities, and oldest upper plates. In these cases, the strong gradients in lithospheric thickness played a role in where lithospheric drips developed [Le Voci et al., 2014].
It is well established that SSC from a 3-D sheared plate preferentially forms longitudinal Richter-rolls (i.e., with their axes aligned perpendicular to the trench) rather than the transverse rolls (with axes aligned parallel to the trench) that can be modeled in 2-D [e.g., Richter, 1973;Wirth and Korenaga, 2012]. Longitudinal rolls have shorter onset times than transverse rolls, and can lead to significant complexity in fully 3-D models of a sheared oceanic plate [e.g., van Hunen et al., 2003;Huang et al., 2003;Ballmer et al., 2011]. A recent study by Wirth and Korenaga [2012] predicts that longitudinal rolls, with wavelengths of 100-200 km and temperature fluctuations of 100-150 K, should occur in the mantle wedge, if viscosities are Շ10 18 Pa s. However, the single mode approximation used in their simulations did not allow for explorations of the full 3-D flow geometry. Honda and coworkers [e.g., Honda and Saito, 2003;Honda and Yoshida, 2005;Honda, 2008;Honda et al., 2010] studied wedge flow patterns in 2-D and 3-D, for cases where only a $100 km wide section of the subarc region was hydrated. For viscosities below Շ5 Á 10 18 Pa s, their models predicted a time-dependent pattern of SSC, where longitudinal rolls interchanged their up and down limbs approximately every 2 Myr. This frequency is comparable to volcanic periodicity in Central Honshu, at arc locations that differ by around 50 km, which is similar to a half-roll wavelength in their models.
In this paper, we extend our systematic 2-D study of SSC in the subduction zone mantle wedge [Le Voci et al., 2014] to 3-D. We explore a range of (dry to hydrous) mantle viscosities, subduction velocities, and upper plate ages. Our 3-D model setup is intentionally simple, with subducting plate motion and the coupling between downgoing and overriding plates prescribed kinematically. We also focus on end-member cases of uniform wedge hydration, making our setup directly comparable to the models of our 2-D study, which allows us to focus on the impact of the additional dimension (section 2). We analyze the resulting 3-D Geochemistry, Geophysics, Geosystems 10.1002/2015GC006125 SSC flow styles and how they vary with different controlling parameters (section 3). Finally, in section 4, we show that 3-D SSC leads to larger spatial and temporal variations in temperature and lithospheric thickness than our 2-D cases, and this more readily produces the conditions required for dehydration and melting.

3-D Model Setup
Our model setup is comparable to our 2-D study [Le Voci et al., 2014], using identical physical parameters that are summarized below. Prescribed slab velocities drive fully-dynamic flow in the mantle wedge beneath an upper plate that is free to destabilize. We solve the Stokes and energy equations, assuming an incompressible, Boussinesq formulation, using the Fluidity computational framework [Davies et al., 2011], which has been carefully validated [e.g., Davies et al., 2011;Kramer et al., 2012;Le Voci et al., 2014] and applied to a range of geodynamical problems [e.g., Hunt et al., 2012;Garel et al., 2014;Davies and Rawlinson, 2014]. The solution strategies employed are identical to those of our 2-D study [Le Voci et al., 2014].

Geometry, Boundary, Initial Conditions, and Material Properties
Our model setup is illustrated in Figure 1. We examine a set of cases in a computational domain that is 400 km deep (z direction) and extends 700 and 1400 km in the trench-perpendicular (x) and trench-parallel (y) directions, respectively. The 50 km thick subducting plate, which is always 50 Myr old at the trench, spans the domain's full trench-parallel extent. Mesh spacing varies from a minimum of 1 km in the wedge corner to 5 km at the domain's base.
Horizontal velocities are prescribed in the incoming plate for 50 km prior to subduction. The slab then follows a down-dipping circular arc to a depth of 75 km, below which it dips at a constant 508 angle. The upper plate is not fully fixed or rigid: it is free to evolve self-consistently in response to the local thermal structure and flow-field [e.g., Kelemen et al., 2003], aside from: (i) at its surface, where we impose a no-slip boundary condition; and (ii) in a curved prism-shaped region above the subducting plate, where velocities are fixed to zero to a depth of 80 km, thereby yielding the so-called ''cold nose,'' which is consistent with observational constraints [e.g., Wada and Wang, 2009;Syracuse et al., 2010]. Below this depth, the subducting slab and mantle wedge are fully coupled, in a manner consistent with the D80 model of Syracuse et al. [2010]. Geochemistry, Geophysics, Geosystems

10.1002/2015GC006125
Vertical boundaries at x5250; 650 km are stress-free, except at x 5 650 km to a depth where temperature reaches 99% of mantle temperature, where a free-slip boundary condition is imposed, to preclude in/outflow at lithospheric depths. At the domain's base, an outflow boundary condition, equal to the subduction velocity, is prescribed on the wedge side of the slab, with a stress-free boundary condition imposed for the sub-slab basal surface. Free-slip boundary conditions are specified at the domain's front and back faces (y 5 0, 1400 km).
Thermal boundary conditions comprise a surface temperature, T s , fixed to 273 K, and zero heat-flux conditions at the domain's base, in addition to the domain's front and back faces. On the x5250; 650 km sides, temperatures are fixed to an error function: where d is the depth (4310 5 2z), T 0 is the reference mantle temperature (1623 K), t plate is either the subducting slab age (t slab ) or the upper plate age (t upper ), while j is thermal diffusivity. Equation (1) is also utilized in defining initial temperature conditions, where, for all cases, t plate 5t slab on the downgoing plate side of the trench, and t plate 5t upper on the upper plate side of the trench.
Material properties are identical to those of Le Voci et al. [2014], with key parameters given in Table 1. Standard values are used for equation of state parameters and these do not vary spatially. Consistent with our current understanding of shallow-mantle rheology [e.g., Karato and Wu, 1993], a composite rheology is utilized, with diffusion (diff) and dislocation (disl) creep viscosities given by: (2) where the subscript X is D for cases with a dry rheology and H for cases with a hydrated rheology. The exponent n in the power law relationship between viscosity, g disl , and the second-invariant of the strain-rate tensor, _ , defines the stress dependence under dislocation creep. T and P represent absolute temperature (for which an adiabatic gradient of 0.5 K/km is added to our Boussinesq potential temperature solution) and Here, C OH and r represent the water content and water content exponent, respectively [Hirth and Kohlstedt, 2003]. A composite viscosity is obtained by combining diffusion and dislocation creep viscosities via a harmonic mean, with viscosity subsequently truncated at a fixed maximum, g max 51310 24 Pa s (no minimum viscosity is imposed). Subslab viscosities are set to a constant value of 1310 23 Pa s.

Simulations Examined
To unravel the dominant controls on the mantle wedge's 3-D flow-regime and thermal structure, we systematically vary a range of subduction parameters (Table 2). To examine the role of viscosity, we vary wedge hydration between dry, ''damp'' [C OH 51000 H/10 6 Si-representative of subridge mantle: Hirth and Kohlstedt, 1996] and ''very-wet'' [C OH 55000 H/10 6 Si-as an end-member subduction zone hydration case: Karato, 2003;Katz et al., 2003], with hydration assumed to be uniform throughout the wedge in all but one case. ''Wet'' cases, with C OH 53000 H/10 6 Si, displayed similar behavior to the very-wet models. We also vary the subduction velocity, over a representative range (2 cm/yr: slow; 5 cm/yr: intermediate; and 10 cm/yr: fast) [e.g., Lallemand et al., 2005;Seton et al., 2012]. For the upper plate, we consider relatively young (50 Myr) and old (120 Myr) cases, with initial lithospheric thicknesses, as defined by the depth of the 1400 K isotherm, of 67.5 and 105 km, respectively. All simulations are run for 50 Myr, and we examine the spatial and temporal evolution of each, rather than concentrating on steady state solutions.

3-D SSC Styles and Controls
We first focus on the observed styles of 3-D wedge flow, which show significant differences from those predicted in 2-D (Section 3.1). We find that the morphology, wavelength, and temporal evolution of instabilities varies across the parameter space examined. Accordingly, we next quantify how the wedge's flow regime and thermal structure varies with: the level of wedge hydration (section 3.2); subducting plate velocity (section 3.3); upper plate age (section 3.4); and the hydrated region's extent (section 3.5). The simulations examined and the resultant flow regimes are summarized in Table 2.

Richter-Roll Style
We illustrate the main characteristics of 3-D SSC using very-wet 120 Myr old upper plate cases. Figure 2 shows a case with a subduction velocity of 10 cm/yr and illustrates how columnar drips align and coalesce into downwelling ridge-like structures to form longitudinal Richter-rolls [Richter, 1973]. The axes of these rolls are principally aligned perpendicular to the trench, which differs to the transverse rolls (with axes Geochemistry, Geophysics, Geosystems 10.1002/2015GC006125 aligned parallel to the trench) that form in 2-D, where such an alignment is not possible. Transient lithospheric drips, which sometimes extend vertically into the wedge core, can be seen propagating along these ridges, toward the wedge corner, due to background corner flow.
Richter-rolls exhibit distinct characteristics beneath the ''arc'' (i.e., the region above the subducting slab at a distance of % 175-275 km from the trench) and ''back-arc'' (i.e., distances տ 275 km from the trench) regions, as is illustrated in Figure 3, for a case with a subduction velocity of 5 cm/yr. Cross sections at 100 km depth (Figures 3a-3c) illustrate that below the arc region, a set of rolls exhibiting alternating high and low temperatures (variations of 75-150 K; Figure 3a) and positive and negative vertical and trenchparallel velocities (Figure 3b/3c) is formed, at a variable wavelength of $100-150 km.
The morphology of subback-arc instabilities differs to those beneath the arc region, as illustrated by cross sections at 150 km depth (Figures 3d-3f): elongated low-temperature ridges are observed, which are 100-150 K cooler than surrounding material, with clear drips along some ridges (as evidenced by regions of lower temperature along ridges in Figure 3d). Ridge spacing ranges from 150-400 km, while they extend in the trench-perpendicular direction for 200-300 km, occasionally to the edge of the domain.  Figure 1 for size). In Figure 2a, which is viewed from below the mantle wedge, the 1550 K temperature isosurface, colored by vertical velocity, illustrates a series of cold ''ridges'' that are principally aligned in the trench-perpendicular direction. Transient ''drips,'' which sometimes extend into the wedge-core, can be seen propagating along these ridges, toward the wedge-corner. Ridges mark the downwelling limbs of longitudinal Richter-rolls, the influence of which can also be seen at the slab surface, where the 1550 K isosurface protrudes further into the wedge; (b) the 1550 K isosurface, viewed from above, which is colored by trench-parallel velocity, providing an alternative illustration of the longitudinal Richter-rolls. In the example shown, ridges are strongly aligned and extend from the wedge-corner to the edge of the domain, often branching or merging; (c/d) an illustration of the mantle wedge's flow regime, from two different directions-images include the 1550 K isosurface (colored in grey) and stream-tracers, which are colored by (c) vertical velocity and (d) trench-parallel velocity. The example shown has a high subduction velocity of 10 cm/yr and, hence, the stream-tracers shown a clear corner-flow pattern (c). Nonetheless, due to the Richter-rolls and transient drips, material can descend into the wedge-core, without reaching the wedge-corner. Furthermore, although corner-flow persists, there is a significant trench-parallel flow component, with material pulled toward the downwelling ridges from both directions, as illustrated in Figures 2b and 2d.

10.1002/2015GC006125
Subback-arc instabilities partially modulate the location of subarc instabilities, imparting a longer wavelength on subarc rolls and locally enhancing subarc downwellings. This interaction is further illustrated in Figure 4, where temporal snapshots of the case illustrated in Figure 3 are displayed at 5 Myr intervals. The location and expression of instabilities, both beneath the arc and the back arc, are strongly timedependent. The subarc system generally shows more temporal variability, with individual rolls merging and splitting on a 1-5 Myr time scale. As a consequence, the number of subarc Richter-rolls varies from 22 to 28 over time (Figures 4a-4d). Beneath the back arc, instabilities coalesce into larger-scale ridges, mostly aligned perpendicular or subperpendicular to the trench (Figures 4e-4h). These ridges often branch or merge, over 5-10 Myr time scales, while also migrating in a trench-parallel direction. This behavior is similar to what was found in 3-D sheared oceanic plate models by van Hunen et al. [2003].

Influence of Viscosity
We find that viscosity exerts the main control on whether or not SSC occurs, which is consistent with previous wedge flow studies [e.g., Honda and Yoshida, 2005;Wirth and Korenaga, 2012;Le Voci et al., 2014]. Figure  5 illustrates the thermal structure and flow regime for a damp case with a 120 Myr old upper plate and a subduction velocity of 5 cm/yr, at t 5 32 and t 5 38 Myr. In a result that is consistent with our 2-D models, we find that damp cases dominantly exhibit instabilities in the form of cold ripples rather than drips [Le Voci et al., 2014]. Interestingly, in 3-D, these ripples can align into 3-D transverse rolls with a weak longitudinal component, which form beneath the back-arc region and rapidly migrate toward the wedge corner. At t 5 32 Myr, the weak longitudinal rolls are evidenced by minor trench-parallel variations in V y in Figure 5f, with 10-12 rolls observed across the trench-parallel extent of the domain, at a characteristic wavelength of 200-250 km (i.e., larger than in the very-wet cases illustrated in Figure 3). At t 5 38 Myr, the transverse roll propagates into the lower-viscosity subarc region and evolves into a set of dominantly longitudinal Richterrolls, with a clear expression in both the temperature and velocity fields, at a slightly smaller wavelength than below the back arc (Figures 5g-5l). Both arc and back-arc rolls are transient features. Note that, in dry cases, no rolls develop in either the arc or back-arc region, over the entire evolution time of our models.
This variation in flow style with viscosity is best illustrated via temporal snapshots of velocity components (scaled with subduction velocity), as a function of distance from the trench. Velocities at 100 km depth ( The dry case exhibits a clear corner-flow velocity pattern, with only minor trench-parallel flow and variability. At 100 km depth (Figures 6g and 6i), trench-ward (negative V x ) and positive vertical velocities below the arc correspond to flow into the pinch zone at the wedge corner, with flow velocities essentially zero elsewhere (this depth is close to the base of the upper plate lithosphere). At 150 km depth (Figures 7g and 7i), cornerflow corresponds to almost purely trench-ward flow (negative V x , V z % 0) below the back-arc region. At both depths, the slab exhibits similar magnitude positive V x and negative V z flow components.
Damp and very-wet cases display, on average, the same corner-flow velocity patterns, at 100 as well as at 150 km depth (Figures 6 and 7a, 7d). However, temporal and spatial variability increases substantially with increasing levels of wedge hydration, due to the presence of subarc Richter-rolls and transient sublithospheric ripples (Figures 6 and 7b, 7e). Below the back arc (Figures 7b and 7e), vertical velocity variations are largest, indicating drip-like instabilities. While the average velocity patterns and magnitudes from 3-D cases are very similar to the corresponding 2-D cases, the variability, particularly in the very-wet case, is greater in 3-D. For the very-wet case (Figures 6 and 7b), velocity variations exceed slab velocities, which is consistent with the predictions of Wirth and Korenaga [2012].
The distinct characteristics of subarc and subback-arc systems arise due to differences in viscosity between the wedge-corner and below the upper plate. The overarching control on wedge viscosity is the level of hydration. However, as illustrated in Figure 8, the wedge-corner's thermal structure, where isotherms are compressed against the slab's surface and a hot tongue of mantle material ascends upwards into the overriding plate, leads to strong variations in viscosity between the arc and back-arc regions, with subarc viscosities further reduced through high strain-rates and lower pressure/depth (through the activation volume).
Geochemistry, Geophysics, Geosystems Geochemistry, Geophysics, Geosystems 10.1002/2015GC006125 As a consequence, spatial variations in viscosity across the warm mantle wedge can exceed an order of magnitude. Subarc and subback-arc viscosities are shown in Figure 9, where symbols and error-bars denote trench-parallel averages and trench-parallel variability, respectively. Viscosities are presented for cases with a 120 Myr old upper plate (i.e., discussed here and in section 3.3) and for a set of cases with a 50 Myr old upper plate (discussed in section 3.4). To allow for direct comparison between cases, viscosity values are consistently extracted at t 5 45 Myr. This time-frame, however, does not capture the transient rolls that occur in the damp case ( Figure 5) and, hence, no trench-parallel variability is apparent in the viscosity estimates of the damp cases.
We find that viscosities are comparable to those observed in 2-D. Viscosities at the base of the back-arc upper plate (red circles) can be up to an order of magnitude higher than below the arc (blue squares). The observation that, without exception, viscosities are lower below the arc than in the subback-arc region explains: (i) the difference in instability morphology and wavelength in both regions; and (ii) why the subarc system is generally more time-dependent than the subback-arc system.
We note that in our 2-D models, SSC did not occur if mantle viscosities exceeded $5 Á 10 18 Pa s. Although sublithospheric ripples were observed in cases at viscosities of $1 Á 10 19 Pa, these did not detach from the lithosphere's base and had a negligible influence on the wedge's flow regime. In 3-D, however, for cases with thicker upper plates, these ripples form transverse rolls, with a weak longitudinal component, at Geochemistry, Geophysics, Geosystems 10.1002/2015GC006125 mantle viscosities of Շ10 19 Pa s. Previous results, from models that used a single-mode 3-D approach [Wirth and Korenaga, 2012] or limited hydration to the wedge-corner [e.g., Honda, 2008], predicted SSC cut-off values of 1-5 Á 10 18 Pa s. The higher SSC cut-off viscosities predicted herein is likely due to the fact that our models are fully 3-D, while instabilities are able to develop over a larger region than the hydrated wedgecorner defined by Honda and coworkers [e.g., Honda, 2008] (see section 3.5 for further discussion).

Influence of Subduction Velocity
Horizontal cross sections highlighting subarc and subback-arc instabilities for slower (2 cm/yr) and faster (10 cm/yr) subduction velocity cases are presented in Figure 10. As in our 2-D models [Le Voci et al., 2014], we find that subduction velocity does not control whether or not SSC occurs. However, it does affect the style of SSC. In general, increased subduction velocities lead to more prominent back-arc ridges, with a stronger trench-perpendicular alignment, when compared to cases with a smaller subduction velocity ( Figure 10). The trench-perpendicular extent of subback-arc ridges also increases with increased subduction velocity; ridges extend from the wedge corner to the boundary of the domain only in cases where V slab 510 cm=yr. Furthermore, these cases are temporally more stable, with subback-arc ridges migrating and interacting less.
Contrary to theoretical prediction [e.g., Richter, 1973;Wirth and Korenaga, 2012], however, we do not observe a monotonic increase in the temperature anomalies associated with Richter-rolls with increasing subduction velocity. We find that thermal anomalies are most pronounced for slow and intermediate subduction velocity cases (cf. the 1500 K isotherm at 150 km depth in Figures 3d and 10b/10d), which can be understood from our 2-D results [Le Voci et al., 2014]. In slow subduction velocity cases, Rayleigh-Taylor drips have sufficient time to develop, as background mantle flow is insufficient to align and coalesce these drips into elongated ridges. Conversely, at higher subduction velocities, strong background corner-flow shears drips into sheets before they can fully develop. Accordingly, competition between the time available For all subduction velocity cases, subarc rolls form with similar wavelengths of 100-150 km, which is consistent with predictions from the single-mode approximation of Wirth and Korenaga [2012] (note that these wavelengths are measured from the V y field, rather than the temperature field, as this better highlights individual rolls). However, we find that with increasing subduction velocity, the modulating influence of subbackarc ridges on the sub-arc system becomes more prominent, superimposing a second, longer wavelength of 150-400 km on the subarc system (as illustrated in the thermal field of Figure 10c). We note that the subbackarc instability wavelength becomes increasingly dominant below the arc region as simulations evolve.
Further insight into the effect of subducting-slab velocity on the flow regime beneath the arc can be gained by analysing mantle velocities at 100 km depth (supporting information Figure S1). We observe a decrease in the strength of trench-perpendicular and vertical velocity components, with respect to slab velocity, as subduction velocity increases, while transients, which are the expression of SSC, become more prominent with decreasing subduction velocity (cf. supporting information Figures S1b and S1e). At 150 km depth (cf. supporting information Figure S2), with increasing subduction velocity, vertical velocity variations below the back arc become more uniform with distance from the trench, indicating that drips are increasingly sheared into ridges, which is consistent with the cross sections shown in Figure 10.

Influence of Upper Plate Age
Thicker (i.e., older) plates are well known to be more unstable than thinner (i.e., younger) plates [e.g., Davaille and Jaupart, 1994]. We examine cases with a 50 Myr old (i.e., younger and thinner) upper plate and V slab 55 cm/yr, at different levels of wedge hydration. Cross sections are analyzed at shallower depths (80 and 130 km) to best capture the Richter-roll systems below these thinner plates.
In the resulting thermal structure for a very-wet case (Figure 11), prominent subarc and subback-arc SSC is visible. The characteristics of SSC beneath the arc and back-arc regions, however, are not as distinct as that for a 120 Myr old upper plate case. The reason being that wedge-corner erosion by mantle flow is reduced, thus producing a weaker gradient in thickness (and viscosity) between the arc and back-arc regions. Nonetheless, as illustrated in Figure 11, the wavelength of subarc instabilities remains smaller than those beneath the back arc. We note that the amount of wedge-corner erosion is limited by the imposed 80 km decoupling depth. Should the decoupling depth vary with temperature, more substantial differential thinning may occur between the arc and back-arc regions [e.g., Arcay et al., 2006Arcay et al., , 2007, which would enhance the difference between arc and back-arc systems for younger upper plates.  Figure 8b. Note that the viscosity scale is truncated, to highlight the variability within the mantle wedge. The wedge-corner's thermal structure, where isotherms are compressed against the slab's surface and a narrow hot tongue of mantle material ascends upward into the overriding plate, leads to dramatic variations in viscosity between the arc (i.e., the region above the subducting slab at a distance of % 175-275 km from the trench) and back arc (i.e., distances տ 275 km from the trench) regions, with subarc viscosities further reduced through high strain-rates and lower pressure/depth (through the activation volume).

Influence of Hydrated Region Geometry
Studies into the stability fields of hydrous minerals indicate that slab dehydration is likely limited to a depth of $200 km [e.g., Schmidt and Poli, 1998;Hacker et al., 2003]. Although the exact mechanism by which water is transported through the wedge remains uncertain, several of the proposed mechanisms lead to quasivertical migration from the point of release [e.g., Gerya et al., 2006;Zhu et al., 2009;Wilson et al., 2014], suggesting that significant hydration maybe limited to the wedge-corner.
As discussed in our 2-D study [Le Voci et al., 2014], localized wedge-corner hydration introduces a step change in wedge and upper plate strength. This facilitates instability nucleation, but also limits the region over which instabilities can form: instabilities only form if the hydrated corner extends over a distance similar to or larger than the instability wavelength (i.e., տ150 km from the decoupling point in our 2-D models).
To illustrate the effect of such localized hydration in 3-D, we have examined a case with a 120 Myr old upper plate and V slab 55 cm/yr, where a very-wet hydrated corner extends 200 km from the decoupling point, with the remainder of the mantle wedge damp, similar to background mantle that is sampled below mid-ocean ridges [e.g., Hirth and Kohlstedt, 1996]. Figure 12 and supporting information Figure S3 show temporal snapshots of the thermal structure and flow regime for this case.
As illustrated in Figure 12, the subarc region initially exhibits longitudinal Richter-rolls with a constant wavelength of $120-150 km. Such a constant instability wavelength differs to that of the uniformlyhydrated simulation (see Figure 3), where Richter-roll spacing beneath the arc region is modulated by longer-wavelength subback-arc instabilities, which are more subdued in our variably hydrated case.
Similarly to the uniformly hydrated damp case, a sublithospheric transverse roll can be seen forming beneath the back-arc region around x 5 450 km (Figures 12d-12f) with weak longitudinal rolls, of wavelength $200-250 km, developing along its length. As in the uniformly hydrated case, the transverse roll extends across the domain's entire trench-parallel extent, and rapidly migrates towards the wedgecorner, where it sharply enhances the strength of subarc longitudinal Richter-rolls, and also imparts the longer wavelength of its longitudinal rolls on to the subarc instabilities (Figures 12g-12i). The wavelength Geochemistry, Geophysics, Geosystems 10.1002/2015GC006125 of subarc longitudinal Richter-rolls subsequently returns to $120-150 km and their strength decreases somewhat (supporting information Figure S3), although they do not wane to their previous vigor. These trends imply that if only part of the mantle wedge is significantly hydrated, subarc SSC will be strongly time-dependent, with significant temporal and spatial variations in Richter-roll wavelength, morphology, and vigor.

Consequences of 3-D SSC
In 2-D, we found that under hydrated mantle conditions, SSC can thin subarc lithosphere by a few km and subback-arc lithosphere by 10-15 km [Le Voci et al., 2014]. Such lithospheric thinning extended the region where wet melting was possible, from absent (for dry cases) to spanning most of the arc and back-arc region, but disrupted by drips, for hydrated cases. These drips also modified slab-surface temperatures and shifted dehydration boundaries by up to 100 K and 20 km, respectively. In this section, we analyze each of these consequences for our uniformly hydrated 3-D models.

Upper Plate Lithospheric Thickness
We use the depth of the 1400 K isotherm as a proxy for lithospheric thickness. Figure 13 illustrates that, consistent with our 2-D results, lithospheric erosion below the arc and back-arc regions is generally most efficient: (i) in cases with higher levels of wedge hydration; and (ii) for cases with older/thicker upper plates. Figure 10. Temperature cross sections at (a,c) 100 and (b,d) 150 km depth, respectively, for very-wet (C OH 55000 H/10 6 Si) cases with a 120 Myr old upper plate and a subduction velocity of: (a,b) 2 cm/yr; (c-d) 10 cm/yr. All slices shown are at t 5 45 Myr. The corresponding 5 cm/yr subduction velocity case is shown in Figure 3.

10.1002/2015GC006125
For subarc lithospheric thicknesses, a minor, but systematic, difference is observed between the 2-D and 3-D cases (Figure 13a/13c). For dry and damp cases, lithospheric thicknesses below the back-arc region from 2-D simulations are only slightly larger than those from 3-D models (Figure 13b/13d). However, in very-wet cases, the differences between 2-D predictions and the average of 3-D models can be up to 8 km. For very-wet These have a negligible effect on estimates of surface heat flow, but strongly affect melting conditions (see section 4.2 below). We note that trench-parallel variations in lithospheric thickness of up to 7 km are also predicted for damp cases, when transverse rolls enter the wedge-corner and spawn longitudinal Richter-rolls beneath the arc region (these, however, are not visible in the temporal snapshots shown in Figure 13).
The increased efficiency of lithospheric erosion in 3-D is a consequence of several combined factors: (i) as noted above, 3-D SSC is observed at higher viscosities than comparable 2-D cases, implying that 3-D wedge models are more susceptible to SSC, which is consistent with the predictions of Richter [1973]; (ii) the onset time of SSC is marginally reduced in 3-D, in the presence of sublithospheric shear, which is consistent with the results of Huang et al. [2003]; and (iii) the longitudinal Richter-rolls observed in 3-D are more vigorous than the transverse rolls observed in 2-D (cf. maximum velocities in Figure 10b/10c). When combined, these lead to more efficient lithospheric thinning in 3-D.
However, even with the additional thinning in the 3-D models, our back-arc lithospheric thicknesses remain greater than the % 60 km back-arc thicknesses inferred from heat-flow and seismic velocities by Currie and Hyndman [2006]. Currie et al. [2008] attributed this strong thinning to additional hydrous plate weakening through viscoplastic mechanisms. Arcay et al. [2006] suggest an alternative mechanism: the development of an intraplate decoupling level associated with the formation of a mechanically weak mineralogical layer that is enriched in water resulting from lithosphere hydration.

Wedge Thermal Structure and Melting
In Figure 14, we illustrate thermal conditions below the arc by plotting trench-parallel temperature averages and ranges, at 80 and 100 km depth, for 50 and 120 Myr old upper plate cases, respectively, alongside wet and damp solidi from Katz et al. [2003] (an adiabatic gradient of 0.5 K/km has been added to our Boussinesq potential temperature solution, which is also the case for Figure 15). Note that at distances beyond 250-300 km, these depths are just below the base of the thermal lithosphere for very-wet cases (at all subduction velocities, and for old as well as young upper plates), but within the upper plate for damp and dry cases, which is reflected in the $50-75 K higher average temperatures for very-wet than for comparable damp and dry cases.
Due to stronger lithospheric thinning, the 3-D models predict elevated mantle temperatures at shallower depths than in our 2-D models [Le Voci et al., 2014]. As melting temperatures increase faster with depth than adiabatic temperatures, this increases melting potential. Average and maximum temperatures below the arc region are highest, relative to local melting temperatures, for cases with: (i) higher levels of wedge hydration; (ii) faster subduction velocities; and (iii) younger upper plates. As predicted in 2-D, subarc wedge temperatures are sufficient to induce wet melting, locally, in parts of all except the dry models (this comparison assumes that locally, wet conditions may permit melting even if hydration is insufficiently pervasive to affect wedge rheology). However, consistent with our 2-D results, in all cases examined, temperatures are too low to induce dry melting.
For very-wet cases (and at certain temporal stages of the damp cases) trench-parallel temperature variations of over $150 K are observed, which exceed the transient 50-100 K variations predicted in 2-D. This will expand the regions where melt pockets may form below the arc and back arc. Such variations in thermal structure and melting may be responsible for the complex seismic velocity structures imaged, for example, beneath the Japanese volcanic arc [e.g., Tamura et al., 2002] and have also been related to the local spacing in arc volcanism [Honda and Saito, 2003;Honda and Yoshida, 2005].

Slab Surface Temperatures
Finally, we consider the effect of 3-D arc rolls on slab-surface temperatures (SSTs, measured at the top of the kinematically defined slab). Figure 15 illustrates the trench-parallel range of SSTs for the different cases examined, alongside MORB dehydration boundaries, the stability fields for hydrous mantle minerals and water-saturated sediment and mantle solidi [Hacker, 2008;Grove et al., 2012].  and thin black lines mark where basalt would retain 0.5 wt % and 0.1 wt % water [Hacker, 2008].
Geochemistry, Geophysics, Geosystems 10.1002/2015GC006125 As in our 2-D models and the studies of Lee and King [2009] and Syracuse et al. [2010], we find that higher subduction velocities, higher wedge viscosities, and thicker upper plates result in decreased SSTs. Trench-parallel variations in SSTs for the 3-D models are most sensitive to the level of wedge hydration: a maximum range of $60 K is observed for very-wet cases, with the range decreasing to Շ20 K for comparable dry and damp cases (although transient ripples and rolls in the damp cases can lead to trench-parallel variations of տ50 K). Peaks and troughs in SSTs occur at a wavelength roughly corresponding to that of the arc Richter-rolls (as illustrated in Figure 2), implying that 3-D instabilities have a more significant influence on SSTs than those observed in 2-D, where the signature of instabilities was not apparent at the slab surface [Le Voci et al., 2014]. Figure 15 illustrates that even a small difference in temperature between hydrated 2-D and 3-D models, trenchparallel variations for a single model, or variations between individual cases with different subduction velocities, could change the depth range of sediment melting or melting of hydrous mantle by several tens of km, due to the slopes of these melting curves, relative to that of the SSTs. However, the observed trench-parallel temperature variations are insufficient to significantly influence the dehydration of crustal material: in all models, the completion of dehydration for MORB material at the slab's surface (taken as the boundary where water content drops below 0.1 wt %) occurs above 90 km depth. As noted in our 2-D study, in addition to controlling crustal dehydration, SSTs can also affect the dehydration of mantle minerals [e.g., Peacock, 1990Peacock, , 1996van Keken et al., 2011], where they are exposed at the slab's surface (e.g., in oceanic core complexes), while temporal variability in wedge temperatures, immediately adjacent to the slab's surface, may influence where released fluids may be first taken up (and subsequently released) by mantle minerals [e.g., Grove et al., 2012]. The range of predicted model SSTs could lead to lateral variations of up to 20 km in the depth where serpentinite and chlorite break down, which may further contribute to the spatial distribution of arc volcanism [e.g., Wilson et al., 2014] and the complex seismic velocity structures imaged beneath some volcanic arcs [e.g., Tamura et al., 2002].

Conclusions
This study builds on previous 2- , we find that SSC is a common occurrence. However, in our 3-D models, prominent SSC occurs at viscosities below $1 Á 10 19 Pa s; this cut-off is higher than that predicted in our 2-D models ($5 Á 10 18 Pa s) and in previous studies [$1-5 Á 10 18 Pa s: Honda, 2008;Wirth and Korenaga, 2012]. Although the influence of water on mantle rheology remains debated [e.g., Fei et al., 2013;Girard et al., 2013], such viscosities are consistent with asthenospheric viscosities estimated from joint studies of glacial rebound and inferences from plate dynamics [e.g., Iaffaldano and Lambeck, 2014], in addition to observations of seismic anisotropy beneath the Pacific basin [Gaboret et al., 2003]. As in 2-D, the exact form of SSC depends on subduction velocity and wedge viscosity. In the unstable cases with close to critical viscosities ($5 Á 10 18 -1 Á 10 19 Pa s), transient transverse rolls, with axes aligned parallel to the trench, can occur; however, longitudinal Richter-rolls, with axes aligned perpendicular to the trench, are the dominant mode of instability, particularly for a strongly hydrated mantle wedge [Richter, 1973;Wirth and Korenaga, 2012].
These Richter-rolls exhibit distinct characteristics beneath the arc region, where the upper plate is eroded by wedge flow, and the back-arc region. Subarc rolls develop at a wavelength of $100-150 km. In the backarc system, Rayleigh-Taylor drips, spawned from the base of the upper plate, are sheared by background corner flow to form long, linear, cold ridges. These ridges, which form the downwelling limbs of larger-scale longitudinal rolls, extend from the subback-arc region into the subarc region, and have a larger spacing (150-400 km) than subarc instabilities, due to higher back-arc viscosities. Both subarc and subback-arc rolls are time-dependent, migrating, interacting, and coalescing with surrounding instabilities, with back-arc instabilities more strongly modulating the location of subarc instabilities in cases with increased levels of wedge hydration, higher subduction velocities, and thinner upper plates.
It is noteworthy that the system naturally forms instabilities of distinct characteristics in the subarc and subbackarc regions, even without a limitation on the distance to which the wedge is hydrated by the downgoing plate. Limiting the extent of the hydrated region does amplify these morphological differences between the arc and Geochemistry, Geophysics, Geosystems 10.1002/2015GC006125 back-arc systems, and the wedge also shows more substantial transient fluctuations. We note that the development of distinct arc and back-arc systems does depend on the extent of differential upper plate erosion between the wedge corner and the back arc, which, in the case of uniform hydration, is controlled by the decoupling depth between the upper and subducting plates, a depth that was kept fixed in our models, but likely shifts with temperature, hydration, and/or is controlled by segregation of lighter materials into the wedge corner [e.g., Wada and Wang, 2009;Arcay et al., 2006Arcay et al., , 2007Honda et al., 2010;Magni et al., 2014].
Our study demonstrates that 2-D models, in many ways, provide a good approximation to the average of 3-D models. However, lithospheric thinning is marginally enhanced in 3-D, which is sufficient to enlarge the region over which melting can occur. In addition, subarc Richter rolls lead to trench-parallel temperature and lithospheric thickness variations of $150 K and $5220 km, respectively. These fluctuations provide a potential mechanism for explaining trench-parallel variations in heat-flow, seismic structure, and magmatism. The wavelength of subarc Richter-rolls predicted herein is larger than the common spacing between volcanic centers [e.g., de Bremond d' Ars et al., 1995]. Nonetheless, as was proposed by de Bremond d' Ars et al. [1995], smaller spacing could result from strong time dependence in the position of high-temperature regions, as is predicted in our models.
On Earth, wedge thermal structure is likely further affected by a heterogeneous distribution of volatiles and strong gradients or steps in lithospheric thickness. As verified in Le Voci et al. [2014], spawning of sublithospheric instabilities can be facilitated by the presence of such features. In addition, wedge flow patterns and temperatures may be modified by 3-D slab geometry [e.g., Kneller and van Keken, 2008], flow around slab edges [e.g., Kincaid and Griffiths, 2004], and compositional and melt buoyancy [e.g., Gerya et al., 2006;Zhu et al., 2009]. The addition of such complexities, to simulations and analyses like those presented herein, is an important avenue for future research.