The Role of Mantle Residuum in Subduction-Driven Circulation: Implications for the Cascades System

This paper investigates the role of upper mantle heterogeneity on circulation. We utilize a kinematic laboratory subduction model that drives multiple styles of three-dimensional, time dependent mantle circulation that interacts with residuum left behind after a major melt production and extraction event. Models of mantle-residuum interaction highlight the impact of melt-induced mantle chemical heterogeneity can have on return flow in convergent margins. Wedge chemical heterogeneity and subsequent melting influence mass and energy transport to the Earth’s surface. These models are designed to represent aspects of the Cascadia subduction zone of the Pacific Northwest U.S., a geologically challenging area where the ~20 Ma Columbia River/Steens flood basalt (CSFB) event transitions to the age-progressive volcanic tracks along the High Lava Plains (HLP) and Snake River Plain-Yellowstone (SRP). This region is a focal point for debate over the applicability of mantle plume models. A series of lab models have tested plume or non-plume explanations for the CascadesYellowstone volcanism. We report on non-plume results that indicate that a nonplume explanation, in which trench-normal volcanic tracks form naturally from deformation of residuum from a CSFB-type event. Melt generating patterns show that our models readily explain bimodal melting beneath the HLP and SRP regions and the westward progression of melt beneath the HLP. However, replicating the ageprogressive melting patterns under the SRP remains a challenge for these experiments.

location is on the slab centerline (star A). In the Cascadia cases (Table 1)

Introduction
Melting models tend to focus on the spatial and temporal patterns of surface volcanism, while less attention is paid to the residue of melting that continues to evolve in the mantle. The concept of mantle residuum has a long history in the geosciences literature, particularly in terms of continental tectosphere, where extraction of volatiles in partial melt creates a local separation between geotherm and solidus, driving up residuum viscosity (and buoyancy) relative to ambient mantle [Jordan, 1975[Jordan, , 1978Pollack, 1986;James et al., 2001]. While there have been some attempts to model time-varying evolution of mantle residuum in plate-driven flow fields [e.g., Humphreys et al., 2000;Smith et al., 2009], the phenomenon remains poorly understood. Our goal is to characterize the general aspects of flow-residuum interactions in subduction settings, with specific focus on the Cascades subduction system in the Pacific Northwest U.S.
The Cascadia subduction zone is a particularly complicated tectonomagmatic region. Most notable geologic events of the Cascades include the ~20 Ma creation of the large igneous province Columbia River/Steen flood basalts (CSFB) and subsequent (~12 Ma to present day) opposite age-progressive rhyolite tracks: Snake River Plain (SRP) and High Lava Plains (HLP) [Geist and Richards, 1993;Camp and Ross, 2004] (Fig. 1). Sinking of oceanic lithosphere with rollback and a mode of plate steepening, back-arc extension, lithospheric basal topography, and a mantle plume have each been implicated in the development of Northwestern U.S. volcanic patterns over the past 20 Ma [e.g., Cross and Pilger, 1978;Christiansen et al., 2002;Camp and Ross, 2004;Jordan et al., 2004;Hooper et al., 2007;Long et al., 2012;Kincaid et al., 2013]. One explanation for the time-progressive HLP/SRP volcanism invokes a CSFB residuum interacting with upper mantle circulation to account for the surface magmatic trends [e.g., Humphreys et al., 2000;Smith et al., 2009] However, the timing of that volcanism in our simplified model is not easily matched to the age-progressive magmatic patterns of HLP/SRP.

Geological Constraints
The Pacific Northwest has a long tectonic and geologic history [Humphreys and Coblentz, 2007], but our focus is on the last 20 Ma of magmatic evolution. Three volcanic features dominate this time period: the Columbia River/Steens Flood Basalts (CSFB), the High Lava Plains (HLP), and the Snake River Plain (SRP) [Geist and Richards, 1993;Camp and Ross, 2004] (Fig. 1). The CSFB is a voluminous pulse of basaltic melting and output beginning after 20 Ma and persisting beyond 15 Ma [Swanson et al., 1979;Carlson and Hart, 1987;Long et al., 2012]. Flood basalts surfaced through generally north-south trending dike swarms over a broad north-south region extending between ~42°N and ~48°N latitude, and lying roughly 500 km away from the trench and coincident with proposed region of plate extension [Cross and Pilger, 1978;Eaton, 1984;Wells and Heller, 1988]. Rough estimates for magma output during the flood basalt period range from ~200,000 -450,000 km 3 [Tolan et al., 1989;Camp and Ross, 2004]. Remarkably, it is estimated that ~98% of the magma volume was erupted in the first 2 Ma of activity [Waters, 1961;Swanson et al., 1979;Tolan et al., 1989;Hooper et al., 2007].
The challenge for geochemical/geodynamic modeling of the post-20 Ma Pacific Northwest lies with the two, more recent east-west trending volcanic chains (HLP/SRP) lying to the south of the CSFB complex. The SRP volcanic track is less voluminous than the CSFB and records a striking age progression in rhyolitic volcanism [Armstrong et al., 1975;Humphreys et al., 2000;Pierce et al., 2000;Camp and Ross, 2004;Jordan et al., 2004;Smith et al., 2009]. This track is generally described as beginning near the McDermitt complex (~16 Ma) near the Oregon-Nevada border [Malde, 1991] or the south-western corner of Idaho (~12 -14Ma) and extending north-eastward to the present day Yellowstone hotspot [Camp and Ross, 2004], with an estimated average eastward shift in felsic volcanic centers of 3 -4 cm yr -1 [Leeman, 1982]. SRP volumetric rhyolitic output is estimated at 10,000 -30,000 km 3 [Ellis et al., 2013], an order of magnitude less than the CSFB. Important constraints on SRP formation come from seismic studies showing strong, low P-and S-wave velocity signals aligned with the SRP trend that are confined in depth (~200 -300 km) and north-south extent (200 km) [Schutt and Humphreys, 2004;Jordan, 2005;Waite et al., 2006;. Geophysical estimates attribute this confined shallow mantle structure to a modest excess temperature (50 -120°C), 1 -2% partial melt and a minor (0.25 mg s -1 ) buoyancy flux [Waite et al., 2005;Smith et al., 2009].
Similar to the SRP, the HLP track records smaller volcanic output than the CSFB and exhibits an approximate trench-normal age progressive trend in rhyolites [MacLeod et al., 1975;Christiansen and McKee, 1978;Draper, 1991;Jordan et al., 2004]. Small volumes of rhyolitic magma, ~ 1,000 -1,250 km 3 [Ford et al., 2013], have erupted continuously from ~10 Ma to present day. The interesting twist is that the HLP felsic production youngs to the west, or opposite the SRP trend, at rates of between 3 cm yr -1 (11 to 5 Ma) and 1 cm yr -1 (5 Ma to present) [Jordan et al., 2004].
HLP morphology is similar to the SRP with a less than 200 km wide track-normal expression in surface volcanism [Draper, 1991;Jordan et al., 2004;Meigs et al., 2009]. The mantle wedge beneath the HLP exhibits a somewhat confined slow seismic velocity anomaly that is weaker (2 -4% ) and less continuous along track than the SRP [Long et al., 2012;Wagner et al., 2012]. The stronger velocity reductions are resolved to a shallow 100 km thick section that appears laterally discontinuous beneath the HLP and broadens towards the western end of the HLP from 100 km to 200 km.

Melting Models
Though widely studied, the tectonomagmatic history of these three spatialtemporal events remain poorly understood and greatly debated [Humphreys et al., 2000;Hooper et al., 2007]. All modes of mantle melting and crustal growth have been proposed to explain the striking volcanic patterns in the Cascades, including subduction driven melting caused by mechanics of the Juan de Fuca plate beneath the North American plate [Faccenna et al., 2010;Liu and Stegman, 2012], ridge melting due to back-arc extension [Carlson, 1984;Carlson and Hart, 1987;Christiansen et al., 2002;Long et al., 2012], and a Yellowstone mantle plume [Morgan, 1972;Jordan et al., 2004;Smith et al., 2009;Kincaid et al., 2013]. Despite the complexity, this area provides an excellent setting for testing a spectrum of models for the style of upper mantle convection in convergent margins and the expected surface expression of related melt generation/transport processes.
Plume and non-plume models have been developed to reconcile the voluminous CSFB (~20 Ma), the age progressive (12 Ma to present) SRP volcanism that terminates at Yellowstone and the opposite, westward trending HLP volcanic track of central/eastern Oregon. There is presently intense ongoing debate on the plume model with the Pacific Northwest at one of the focal points for this dispute [Foulger, 2003[Foulger, , 2005Anderson and Natland, 2005;Foulger et al., 2005]. A plume, or thermally buoyant mantle upwelling with a flood basalt producing head and long-lived tail for generating an eastward progressive track, is commonly cited as the key dynamic process for CSFB and the age progressive rhyolitic volcanism of SRP [Richards et al., 1989;Camp and Ross, 2004;Pierce and Morgan, 2009;Smith et al., 2009]. But the equivocal HLP, whose rhyolite lavas young opposite North American plate motion, and the north-south offset between CSFB and SRP is inconsistent with the simple plume model [Christiansen et al., 2002].
A non-plume model has also been considered to explain the volcanic patterns of the PNW. Supporters of a non-plume theory invoke plate motions to trigger largescale melting events, like the CSFB [Faccenna et al., 2010;Liu and Stegman, 2012;Long et al., 2012]. Geodynamic 3D models have shown that changes in slab motions (e.g. subduction initiation, rollback initiation or acceleration, and slab steepening) can result in an upward mantle circulation (i.e. a proxy for decompressional mantle melting) in the wedge [Kincaid and Griffiths, 2004;Long et al., 2012]. Long et al. [2012] propose that the voluminous CSFB melting event was a direct consequence of  [Humphreys et al., 2000;Smith et al., 2009]. Results are presented from a series of geodynamic laboratory experiments designed to test the importance of the CSFB mantle residuum on 4D circulation, transport, and melt production.

Apparatus
We examine the evolution and interaction of mantle chemical heterogeneity, or residuum, in 4D subduction-driven flow using a kinematic laboratory model [e.g., MacDougall et al., 2014]. The mantle is modeled using glucose syrup, i.e. a Newtonian viscous fluid, in a transparent Plexiglas tank (150 cm long × 75 cm wide × 45 cm deep). Similar to Hall and Kincaid [2001], the subducting slab is modeled by reinforced continuous rubber belts moving around rollers (Fig. 2a). A large percentage of the wedge superadiabatic thermal gradient is locked up in the subducting plate.
Thus, the downgoing slab is expected to be the dominant driver of forced convection in the upper mantle. While dynamic subduction of a tabular plate is ultimately the goal for wedge circulation studies, there remain a large number of poorly constrained processes in these models (e.g., rate of convection and upwelling in the mantle wedge, shear stress on the base of the lithosphere, etc.). We choose kinematic models because forced convection of the wedge can be generated in a reproducible fashion.
We employ kinematic subduction modeling where precise downdip (U D ) and translational, i.e. trench retreat, (U T ) motions of the slab are prescribed to mimic slab motions observed in 3D, time-evolving dynamic subduction models [Kincaid and Olson, 1987;Griffiths et al., 1995;Funiciello et al., 2003Funiciello et al., , 2006Schellart, 2004]. An advantage to using kinematic subduction modeling is that each experiment is repeatable and controllable. In this way we are able to characterize the wedge response to mantle heterogeneity in 3D wedge flow fields when key parameters in plate mechanics are varied. The glucose syrup used to model both upper mantle and residuum is a Newtonian viscous fluid with temperature dependent viscosity, described by Olson and Kincaid [1991]: (1) where µ and T are dynamic viscosity (Pascal-seconds) and temperature (°C), respectively. The volume of fluid (4.5 × 10 5 cm 3 ) representing the ambient mantle is maintained at a constant temperature, resulting in isothermal density of 1.42 g cm -3 and dynamic viscosity of 134 Pa s at 22°C. The mantle wedge is isothermal until it interacts with the mantle residuum, a thermally altered patch of glucose syrup. The volume of the residuum is 200 cm 3 for all experiments, while temperature (T R ), thereby viscosity (µ R ), is varied. For visualization purposes the mantle heterogeneity is marked with neutrally buoyant red tracers (beads, whiskers) and micro-bubbles that distinguish it from the ambient mantle fluid, which is laced with white whiskers.
These passive Lagrangian flow tracers move with the fluid providing velocity information in both space and time throughout the course of the experiments.
We model three distinct styles of overriding plate in our experiments. In the first setup, we assume a decoupling zone exists between mantle and overriding lithosphere. These experiments have a free slip mantle surface, in which the mantle surface moves independently of upper plate motion. Our second set of models use a strongly coupled mantle-lithosphere boundary to create a no slip mantle surface condition. In these cases a uniform overriding plate (OP) is modeled using a thin, transparent Plexiglas plate that migrates with trench motion. The OP couples with the underlying wedge fluid, imposing a vertical shear flow under the plate (Fig. 3a). In the third set of overriding plate models we simulate extension in the back-arc (BAE) using mylar sheeting along the fluid surface (Fig. 3b). The extension center moves at the prescribed trench migration rate (U T ). As the trench retreats, mylar spools from a reel to the extensional axis in the back-arc. The mylar couples with the underlying fluid east of the extension center and imposes a zero horizontal surface velocity on the wedge surface.

Scaling
The length and time scales of these kinematic models, representing dynamic sinking modes, are scaled to the mantle using a dimensionless ratio of advective to diffusive heat transport, written as the Péclet number, where length scale, L, is defined as the width of the trench (L lab = 25 cm, L mantle = 1625 km) and the thermal diffusivity (κ) values for lab and mantle are 10 -3 and 10 -2 cm 2 s -1 , respectively. Time and velocity scales are calculated by relating laboratory and mantle Péclet numbers. Thus, a U D of 1 cm min -1 in the lab relates to 0.8 cm yr -1 when scaled to the Earth's mantle and 1 min in the lab corresponds to ~8 Ma of geologic time.

Experimental Methods
The experimental method for melt residuum experiments has been developed through the running of over 60 distinct cases. Experiments begin after the placement of mantle heterogeneity on the surface of the upper mantle wedge, i.e. our experiments begin after a major melting event resulting in creation of the flood basalt residuum.
A residuum volume of ~200 cm 3 scales to an initial melt residuum volume of 5.5 × 10 7 km 3 . We assume a residuum volume ~10 times larger than the combined estimated erupted and unerupted magmas extracted from the mantle during the CSFB melting event. Rough estimates for magma output during the Cascades flood basalt period range from ~200,000 -450,000 km 3 [Tolan et al., 1989;Camp and Ross, 2004]. We also consider that the solidified magma intrusions beneath a large igneous province can represent a substantial unerupted igneous volume [Crisp, 1984;Walker, 1993;White et al., 2009]. Thus, when determining the initial size of our CSFB melt residuum, we estimate a 10% mantle melt fraction represented by surface lavas and unerupted, underplated magmas.
The laboratory residuum is made from the same glucose syrup as the ambient wedge fluid. Temperature differences are used to represent a wide range in residuum viscosity contrasts. Heterogeneity temperature, i.e. viscosity, is an experimental variable. To create the low or high viscosity heterogeneity we heat or chill the glucose syrup, respectively.
After the heterogeneity is brought to the appropriate temperature, it is placed in the wedge. Location in map view is most often centered in the wedge about the slab centerline and ~14 cm (~900 km) from the trench (e.g. location 'A', Fig. 2c). The heterogeneity is approximately oval shaped, long in the trench-parallel direction, ~20 cm (~1300km), and constrained in the trench-normal direction, ~5 cm (~325 km).
While we can control heterogeneity volume, temperature, and general location, once the residuum is placed in the wedge, natural viscous effects take place prior to experiment initiation. Our low viscosity heterogeneity has a viscous puddle spreading speed of ~0.4 cm min -1 when placed in ambient mantle fluid. So naturally less viscous heterogeneities spread more and tend to occupy more surface area than more viscous heterogeneities when placed in the tank. Likewise, less viscous heterogeneities extend to shallower depths (~2.5 cm, corresponding to ~165 km in the mantle), while more viscous heterogeneities reach to deeper depths (~ 3.5 cm, ~225 km). Thus, there are slight initial morphology differences at the initiation of each experiment.
In the following cases, and throughout the paper, we use the terminology "weak" or "strong" heterogeneity to relate the strength, or viscosity, of the residuum (µ R ) to that of the surrounding mantle (µ M ). The viscosity ratio, , is used to determine whether the residuum is weak (µ * < 1) or strong (µ * > 1). We assume that the strength of the flood basalt residuum is determined by the mantle's ability to retain partial melt and volatiles after the CSFB event. In the case of strong residuum we assume eruption of melt and volatiles to the surface results in a more viscous chemically altered heterogeneity [Jordan, 1975[Jordan, , 1978Pollack, 1986;James et al., 2001]. Following a similar logic, a low viscosity heterogeneity results from the mantle's capacity to stably hold some partial melt at shallow mantle depths [Hirschmann, 2010]. Though P and S wave velocity anomalies can detect the presence of retained melt in today's upper mantle [e.g., Toomey et al., 1998;Villagómez et al., 2014], dV P and dV S data only offer a present day snapshot. Since we are unable to determine the residuum viscosity resulting from a melting event that occurred millions of years ago, e.g. CSFB at ~20 Ma, we model both end members (µ * < 1 and µ * > 1) in each set of experiments.
Assuming that the residuum creation occurs at ~20 Ma and that the experimental time begins at t * = 0, we monitor the wedge flow-residuum interaction through present day (0 Ma Reasonable values for the laboratory glucose syrup are ρ = 1.42 g cm -3 , D = 40 cm is the depth of the tank fluid representing the mantle, u = 8 cm min -1 , and µ = 134 Pa s.
The Reynolds number for our laboratory model is calculated to be Re = 5.6 × 10 -3 .
Thus our analogue mantle fluid with Re < 1 is largely unaffected by inertia and flow is laminar.
With the high viscosity heterogeneity present in the shallow wedge, the system is unstable and convects. Using the formula for the onset time of a residuum-ambient fluid boundary layer instability from Olson [1990]: where Ra c is the critical Rayleigh number, υ is the dynamic viscosity, g is gravity, ΔT is the maximum temperature rise between the ambient fluid and residuum, and κ and α are thermal diffusivity and thermal expansivity of the lab fluid, respectively.
Reasonable values of Ra c = 3 × 10 3 , υ = 10 3 cm 2 s -1 , g = 9.8 m s -2 , ΔT = 9 K, κ = 10 -3 cm 2 s -1 , and α = 4.5 × 10 -4 K -1 , gives τ ~ 7 hours, the amount of time for Rayleigh-Taylor instability onset on the underside of the most viscous residuum. Using the scaling indicated by the Péclet number we determine that the diffusion characteristic timescale for corn syrup with thermal diffusion distance of 1 cm is ~17 minutes, almost twice the length of a typical experiment. Subsequently, residuum buoyancy and thermal diffusivity are second order factors over the duration of an experiment.
Our length scale is defined by the width of the trench (~1625 km) which is larger than the current width of the Juan de Fuca trench, ~1200 km. Plate tectonic reconstructions [e.g., Atwater, 1970;Atwater and Stock, 1998] [Atwater and Stock, 1998], resulting in a subducting slab ~2000 km wide. In our lab model, with approximate representation of a system with complex spatialtemporal characteristics, the subducting plate maintains a fixed width throughout the experiment. In order to account for the diminishing width of the Cascades slab over the last 20 Ma, we model an average slab width over this time period of 1625 km.

Data Collection
Data on deformation and entrainment of mantle heterogeneity in subductioninduced flow are collected throughout the experiment duration using high-resolution digital photography. Two cameras are used to take side view and map view Fluid flow patterns represent material pathlines in a fixed trench reference frame.
Throughout this paper, all velocities and flow directions are given in a fixed trench reference frame. For more information on how reference frame affects observed wedge flow fields see the Supplemental Information.

Experimental Results
Using a subset of 15 experiments (Table 1), we investigate the role of upper mantle heterogeneity in subduction systems. We develop qualitative and quantitative relationships between the subducting plate style, overriding plate structure, wedge flow fields, and melt residuum transport and deformation. In each case we characterize the time dependence and distribution of material transport in the wedge and compare these data to the seismologic results, geochemical data, and post-20 Ma surface tectonomagmatic activity of the PNW.

Time Evolution of Residuum: Rollback Subduction
The first set of experiments model a free slip surface, i.e. we model mantle decoupling from the overriding lithosphere. Results show that wedge return flows are highly 3D and the addition of a rheological heterogeneity disturbs basic flows. Fig. 4 highlights the evolution of a rollback experiment where the slab subducts with downdip (U D = 8 cm min -1 ) and translational rollback (U T = 3 cm min -1 ) motions. All images are shown in a fixed trench reference frame.
The weak heterogeneity (Exp. 6, µ * = 0.3) experiences pure shear in transformation from an initial north-south morphology to an east-west morphology after 2.5 minutes (~20 Ma) ( At t * =1 (~1.25 min into the experiment) the leading edge of the residuum is at the trench and will begin subducting from the shallow mantle. Later stage images of this case ( Fig. 5a -b) reveal the long-term fate of weak heterogeneity in the subduction system. After 10 minutes (~80 Ma or t * = 4) most of the initial residuum volume has subducted (estimated ~95%), but a very small amount (~5%) has not flushed from the shallow wedge. Over this time period the low viscosity feature shows continued trench-parallel thinning; with a trench-parallel width of < 1 cm (~ 50 km), extending along the slab centerline. This narrow band of low viscosity fluid remains in the wedge for the duration of the experiment and has > 80 Ma expected lifetime in the shallow wedge.
Changing residuum µ * affects mass and energy transport within the wedge, on both short and long-term timescales. Fig. 4d -f illustrates the ~20 Ma evolution of a strong (i.e., high viscosity) heterogeneity in a rollback experiment (Exp. 5, µ * = 8). All model parameters remain the same from Exp. 6 except for the viscosity contrast (µ * ).
Similar to the weak heterogeneity, slab entrainment translates the strong heterogeneity towards the trench and U T -induced toroidal flow results in north-south compression towards the slab centerline. Unlike the low viscosity wedge feature, which experiences pure shear as it deforms into an east-west trending morphology, the high viscosity feature resists internal deformation and instead folds along the trench-normal symmetry line ( Fig. 4e -f). The strong residuum is more resistant to the deformation flows coming from the north and south edges at rates of ~3.5 cm min -1 (~0.44U D ).
Unlike the low viscosity heterogeneity which compresses down to ~6 cm (~390 km) in trench-parallel width, the more viscous morphology maintains 55% of it's initial north-south length and at 20 Ma has a trench-parallel length of ~11 cm (~715 km).
Another striking dissimilarity between variable µ * cases is the residuum entrainment style. With strong heterogeneity present in the wedge, the entire feature is entrained towards the trench as one block. This feature translates at ~2 cm min -1 (~0.23U D ) which is a slower rate than the weak residuum. By present day it is still ~4 cm (~260 km) from the trench. At later stages of evolution, efficient entrainment of viscous heterogeneity results in ~95% of the heterogeneity flushing from the shallow mantle wedge over 80 Ma (Fig. 5c -d), and 100% by 90 Ma. So while both strong and weak heterogeneities result in ~95% of the residuum being subducted after 80 Ma, the behavior of the remaining 5% of the residuum typifies the varying reactions of the residuum to the flow fields produced in the wedge by the different tectonic drivers.
Circulation in the form of return flows into the innermost wedge is strongly influenced by the presence of viscous heterogeneity. Average instantaneous velocity and particle direction are summarized in 3 distinct wedge regions: slab centerline, northern edge, and mid-region (Fig. 6). For each zone, instantaneous non-dimensional velocity is defined as , where is the slab convergence rate (8 cm min -1 ) and (5) where , trench-normal velocity, and , trench-parallel velocity, are calculated using x, y, and t data collected from digitizing software GraphClick. Normalized particle flow direction (Φ * ) in a fixed trench reference frame is also calculated from this data. Fluid migration direction, Φ * , ranges from 1 to -1, where Φ * = 1 represents westward (trenchward, trench-normal) velocities, Φ * = 0 for southward (trenchparallel) flows, and Φ * = -1 for apparent eastward (away from slab, trench-normal) directional flows. We calculate instantaneous and Φ * for particles located in the ambient fluid and residuum at three times: experiment initiation (t * = 0 or 20 Ma), at 1.25 minutes (t * = 0.5 or 10 Ma), and at 2.5 minutes (t * = 1 or present day); then average over different regions. We confine our analysis to the northern wedge; our simplified model is symmetric about the slab centerline so we observe mirrored flow patterns about the slab centerline in the southern wedge. To isolate the impact of heterogeneity in rollback subduction, Fig. 6 compares average instantaneous velocity, , and direction, Φ * , in the northern (black), mid (blue), and central (red) regions of the wedge (see Fig. 6a for location boundaries) for weak and strong heterogeneity present in the wedge (Exp. 6 and Exp. 5, respectively).
There are basic similarities in wedge circulation, along with key differences, when a weak versus strong residuum is present. Regardless, flow enters the edge of the wedge in the north to south direction and near the slab centerline transitions to fully east to west, trench-normal, approach velocities ( Fig. 6b). On average, there is a gradual transition from trench-parallel to trench-normal flow through the mid-region of the wedge. Trench parallel flow rates entering the edges of the wedge reach ~0.55U D with a weak residuum and slightly reduced rates, ~0.50U D , when a strong residuum is present (Fig. 6c). Average trench-normal flows in the central wedge (red region) are inversely related to distance from the trench though mean velocities are slower, reaching ~0.45U D when µ * < 1 and ~0.35U D when µ * > 1 (Fig. 6e).
There are, however, a number of important differences in detailed flow structure based on viscosity of the residuum. Moving trenchward (from east to west) the greatest increase in U M (+0.30U D ) is seen in the wedge with weak heterogeneity.
With a high viscosity residuum in the wedge, slower slab centerline velocities supply mass to the trench (Fig. 6e). Activated return flows show up in greater U M increases in the northern edge and mid-region when a strong heterogeneity is present in the wedge ( Fig. 6c, d). Values for U M in the strong residuum case show trenchward increases of +0.15U D and +0.05 U D over the weak residuum case in the northern edge and midregion, respectively (Fig. 6c, d).

Effect of an Overriding Plate
Surface flow coupling to an overriding plate is important for full wedge circulation results. Two end member wedge surface conditions exist: fully coupled and fully decoupled. Experiments in the previous section have a free slip surface condition, i.e. decoupled. In the free slip mantle end member condition we assume that the non-Newtonian rheology of the mantle [Hirth and Kohlstedt, 2003] 9). Initial north-south residuum length L 0 = 20 cm, corresponds to 1300 km when scaled to the mantle. A weak heterogeneity experiences the greatest trench-parallel deformation over a 40 Ma (5 minute) evolution. By t * = 2, low viscosity residuum has changed north-south extent dramatically (~0.1L 0 ) with no OP and moderately (~0.5L 0 ) when the wedge couples to an OP. This is consistent with 3D flow rates in Fig. 7. A strong residuum located in the wedge is more resistant to slab induced deformational flows. After 5 minutes (40 Ma) evolution, a high viscosity feature maintains ~0.3L 0 with no OP and ~0.95L 0 with OP. The decrease in trench-parallel deformation in cases with overriding plate is the result of two factors: (1) introduction of an overriding plate dampens U T -induced toroidal near-surface flows (Fig. 8b, c); and (2) residuum, especially µ * > 1, viscously couples to the base of the OP and resists deformation (Fig.   9b).
Results indicate that slab-induced upper mantle circulation can deform and entrain wedge heterogeneity. We've shown that a weak heterogeneity decoupled from the overriding plate most easily deforms from a trench-parallel (north-south) morphology to a trench-normal (east-west) morphology, over 20 Ma time scales.
However, an additional first order feature of volcanism needs to be considered, namely melt production patterns. These include volumetric trends and age progressions in the lavas, e.g. the CSFB event produced 200,000 -450,000 km 3 [Tolan et al., 1989;Camp and Ross, 2004] while the SRP eastward trending and the HLP westward trending rhyolitic volumetric output is estimated at 10,000 -30,000 km 3 [Ellis et al., 2013] and 1,000 -1,250 km 3 [Ford et al., 2013], respectively. To access vertical velocity, we consider flow and deformation patterns in side view for cases of different µ * and surface condition (Fig. 10).
Side view analysis provides insight on the slab entrainment and transport of µ * residuum from shallow mantle to deep (Fig. 10 a, b). At t * = 1, the less viscous heterogeneity resembles a thin vertical sheet that runs along the slab centerline (e.g., Fig. 4c) and extends from mantle surface to depths beyond the shallow mantle (Fig.  10a). Comparatively, analysis of the vertical plane with strong residuum present shows that the feature is fully entrained and the entire mass approaches the trench (Fig. 10b), similar to the map view observation (e.g., Fig. 4f). To compensate for slower shallow mantle velocities when strong heterogeneity is present (Fig. 6), deeper mantle is sourced for subduction. In our x-y plane velocity analysis (Fig. 6c -e), strong residuum in the wedge results in slower circulation throughout the shallow mantle horizontal plane (z ~1 -2 cm). In the x-z plane, the greatest sinking velocities, up to 0.5 cm min -1 , occur when a high viscosity feature is present in the wedge (Fig. 10b). In comparison, weak heterogeneity experiences greater flow velocities in the shallow wedge x-y plane (Fig. 6c -e) and slower sinking velocity, ~0.3 cm min -1 , in the vertical plane (Fig. 10a). Thus, the relative intensities of entrainment and sinking velocities supplying the mass for subduction varies based on µ * heterogeneity present in the shallow mantle wedge.
Residuum couples to the no slip surface when an overriding plate is introduced (Fig. 10c, d) producing a significant decrease in sinking velocities. Weak heterogeneity tends to remain in the shallow mantle wedge, adhered to surface plate, until coupling occurs with the subducting slab (Fig. 10c). More viscous heterogeneity strongly adheres to the underside of the OP and experiences limited deformation and entrainment over 2.5 minutes or 20 Ma when scaled to the mantle (Fig. 10d). There is slight shear-induced upwelling detected beneath and on the western side of the residuum in each case. Vertical rise velocities up to ~0.2 cm min -1 (weak heterogeneity) and ~0.1 cm min -1 (strong heterogeneity) are observed along the deep and western residuum-ambient fluid interface. However, no evidence of vertical motion exists on the eastern side of the heterogeneity.
Long-term flux of rheological heterogeneity to the deep system is highly dependent on µ * and fluid surface coupling conditions (Fig. 11). With no OP, the estimated volumetric flushing time for weak and strong heterogeneity is approximately the same (Fig. 11a). A strong residuum in a system with a decoupled fluid surface (e.g., Exp. 14), µ * is strongly entrained down with the subducting slab.
98% of the heterogeneity is flushed from the shallow wedge to depths greater than ~3 cm (~200 km) by t * = 3, corresponding to 60 Ma when scaled to the Earth (Fig. 11a).
The remaining 2% is strongly entrained, meaning the high viscosity residuum will completely flush from the wedge within the next 5 Ma resulting in a total shallow wedge lifetime of ~65 Ma. The majority of the low viscosity residuum volume subducts as well (~90% in 60 Ma) (Fig. 11a). However, a thin vertical sheet of residue is not strongly entrained (e.g. Fig. 5a -b). Late stage results show that low viscosity features, though thin, are a long-lived (over 80 Ma) source of chemical heterogeneity in the wedge.
Deformation and entrainment patterns of wedge heterogeneity alter drastically when an overriding lithospheric plate is added to the model (Fig. 11b). Over 60 Ma, strong residuum lying beneath an OP experiences slight deformations and entrainment.
But due to strong viscous coupling with the overlying lithosphere, only 6% of the residuum volume has downwelled below the shallow mantle. ~83% of the weak residuum with fully coupled mantle surface has flushed by t * = 3. The remaining 17% exists in a thin horizontal sheet that is coupled to the no slip surface. Here, and in the case with no OP, low viscosity heterogeneities are capable of long lifetimes in the mantle wedge. Thus weak residuum has the greatest potential to impact wedge chemistry over long timescales.

Effect of Back-Arc Extension
Results so far show how simple plate motions give rise to evolving 3D flow fields that affect residuum deformation and evolution, but so far residuum cases have not produced complex upwelling patterns. Although our decoupled versus coupled wedge surface results show strong differences in vertical velocity, these upwelling patterns do not match a number of key characteristics of Cascades, e.g. HLP/SRP thin tracks with opposite age progressive lavas. An essential condition of the Cascadia subduction system that we will be modeling is the introduction of back-arc extension (BAE) above the mantle wedge [Wernicke et al., 1988]. The presence of an extensional axis provides a dramatic return flow regime boundary for material reaching the areas beneath the arc volcanics and the wedge apex. The subducting slab does not entrain residuum and ambient fluid located east of the extension center as both trench and extensional axis translate westward at rates of 3 cm min -1 (Fig. 12). Heterogeneities adhere to the OP, which is stalled relative to the system rollback, and are left behind at ~ 2 cm min -1 (Fig. 12b, e). Wedge material located beneath the AP, i.e. trenchward of the BAE axis, is entrained towards subduction. East-west extension (in the trench-normal direction) occurs along the length of the heterogeneity due to the increasing gradient in trench-normal approach rates.
Horizontal velocities for a weak heterogeneity located beneath the AP reach ~3 cm min -1 . Between 10 and 0 Ma (t * = 0.5 -1), the low viscosity feature separates into two unconnected entities (Fig. 12c). The space between the two discrete features grows with time as the subduction system moves west and the detached residuum beneath the OP is left in the east. Exp. 18 illustrates the strong heterogeneity stress guide effect. Due to viscous coupling within the strong heterogeneity, slab-ward advection of residuum beneath the AP is limited (Fig. 12d -f). Approach velocities of the western side are ~1 cm min -1 , or only ~30% of the rates for the weak heterogeneity case (Fig. 12b). After 20 Ma, the strong heterogeneity remains intact and robust from west to east (Fig. 12f).
Complex 4D flow fields emerge when upper plate extension is modeled. In respect to the trench, wedge material beneath the overriding plate has an average southeast to east flow direction across all zones (Fig. 13b, c). As in previous cases, despite the presence of BAE circulation in the wedge is sensitive to a mantle heterogeneity located beneath the extension location. On the northern edge region (Fig. 13a, black zone) and mid-region (Fig. 13a, blue zone) there is limited residuum presence during the course of the experiments. Between cases average U M and Φ * are similar within each zone. As previously discussed, fast trenchparallel flows exist on the northern edge. Moderate (~0.35U D ), westward trenchoblique particle motion exists on the mid-region (Fig. 13e). There is a greater residuum presence along the slab centerline. In the red zone flows readily reorient to a trench-normal direction from east to west. Wedge flow in this region reorients at a slower rate when strong residuum than when weak residuum is present (Fig. 13c). Slab centerline U M is the same for both cases beneath the OP. Beneath the AP, weak heterogeneity experiences faster average slab-ward velocities (~0.35U D ) than the case with the strong heterogeneity (~0.25U D ) (Fig. 13f).
Investigation of heterogeneity evolution in side view further highlights the importance of µ * features in 4D upper plate-driven wedge (Fig. 14a, b). As we saw in map view, the weak residuum splits into two distinct features while the strong heterogeneity remains intact. Exp. 27 (low viscosity residuum) is distinct in that the ambient fluid-residuum interface is slanted and shallows with time (Fig. 14a). Along the slab centerline, upwelling is detected in the shallow wedge along the residuum interface boundary with the fastest rise velocities (up to ~0.05 cm min -1 ) seen, as expected, beneath the extensional axis (Fig. 14c). The strong heterogeneity is robust directly beneath the BAE location (Fig. 14b). Maximum vertical velocities (up to ~0.05 cm min -1 ) are resolved within and up the eastern side of the viscous plug ( Fig.   14c). In both cases shallow wedge velocities are beginning to downwell (v z < 0) near the trench due to the downward entrainment of the slab. However, vertical rise velocities still exist in the shallowest portion of our domain-averaging window (Fig.   14c, representing the region between mean and one standard deviation above the mean).

Results reveal interesting horizontal versus vertical flow dominant regimes
indicating that upwelling patterns depend not only on residuum presence but also on location within the wedge. It is only along the slab centerline that trench-normal approach velocities, that are less than downdip subduction rate, must be compensated for by an advected vertical velocity. Here, material must be brought up from deep, against gravity, to fill in the mass deficit created by the extension center on the surface. Away from the central wedge flow is predominantly horizontal toroidal flow, with little vertical component. Thus, model predictions indicate that melting will be constrained on trench-normal tracks on the slab centerline.
Analysis of µ * residuum deformation beneath the BAE axis reveals that it is possible to produce age-progressive mantle melting. As previously discussed, we confine the greatest upwelling to central wedge, along the trench-normal, slab centerline track. Relative strength and chemical composition of the residuum beneath the back-arc extension axis will effect upwelling related decompression melting trends (Fig. 14). The interface represents a boundary between the more refractory residuum and ambient mantle. As it shallows and migrates, this boundary controls both the chemical composition and timing of the resulting partial melts (Fig. 15).
Different wedge upwelling patterns occur when heterogeneity viscosity is varied (Fig. 14). Shallowing of residuum and ambient fluid throughout the top 3 cm (200 km) creates a decompression melting favorable environment. An essential point however is that not all wedge upwelling will produce decompression melting. Melt production will be delayed in rising residuum parcels that have experienced prior melt extraction and are more refractory, with a solidus that is shifted towards higher temperatures. Focusing in on the residuum-ambient fluid interface, this results in a different melt onset times for neighboring parcels. Melting will depend not only on a particle's ability to rise but also on its composition, location in P-T space, and wedge potential temperature (Fig. 15a). In our analysis we assume a fertile mantle lherzolite composition for the ambient mantle and a refractory harzburgite composition for the depleted mantle residuum. Using these assumptions, ambient mantle will start to melt at depths between ~75 and 20 km for wedge potential temperatures of 1500 and 1300°C, respectively (Fig. 14d). Refractory mantle heterogeneity will only melt at mantle depths shallower than 20 km in the hottest wedge (T P 1400 -1500°C).
Relative strength of the residuum plays a key role in the wedge spatialtemporal patterns of upwelling. Low viscosity heterogeneity interface shallowing with time can produce an age progressive melting track (e.g., Fig. 15b). Rise velocities beneath a weak mantle heterogeneity decrease with distance from the extension axis Rise of mantle material and shallowing of residuum interface decreases drastically 21 -22 cm (~1400 km) from the slab (Fig. 14a). In the long term, as the extension axis migrates westward, further away from the eastern heterogeneity, upwelling velocities beneath the interface will decrease to zero unless an additional source of upwelling is introduced.
Similar to the weak heterogeneity case, predictions on west and east age progression patterns in melting with a strong heterogeneity present in the wedge are linked to residuum-mantle interface rise in P-T space. The high viscosity residuum (Exp. 18) redirects BAE-induced upwelling around the east and west sides (Fig. 14c).
As the extension axis migrates westward at 3 cm min -1 , vertical flows preferentially rise up the right side of the feature. Consequently, rising of residuum and ambient fluid is strongest in the east. Weaker upwelling is detected west of the BAE axis, as flow near the trench compete with the downward entrainment of the slab. 1 cm min -1 (~0.8 cm yr -1 scaled to mantle) westward progression of residuum is occurring to the west of the extensional axis. However, the eastern side is coupled to an overriding plate that has zero horizontal velocity. Upwelling exists to the east of the BAE location but without an additional driver allowing eastward progression of interface in an absolute reference frame, there appears to be no eastward age progression in this wedge. Only in a fixed trench/fixed BAE reference frame does the eastern edge have relative eastward progression.

Effect of Residuum Initial Position
A north-south (trench-parallel) offset is observed in cases where residuum is initially located north of the slab centerline (Exp. 23 and 30). Trench-parallel translation is a consequence of deformation induced by trench rollback (Fig. 16). Both weak and strong residua are initially located north of the slab centerline (Fig. 2c, location 'B'). Weak heterogeneity migrates south at a maximum rate of ~5.5 cm min -1 .
From 20 Ma to 10 Ma, the residuum's northern most edge has moved ~7 cm to the south, corresponding to ~450 km in the mantle (Fig. 16b). As we saw in previous velocity results, particle trench-parallel velocities decrease towards the slab centerline (e.g., Fig. 13). The southward migration rate of the transported residuum from 10 -0 Ma lessens and the total trench-parallel offset is 550 km after 20 Ma evolution (Fig.   16c). Similarly to previous cases with BAE, northern wedge residuum morphology experiences trench-normal stretching over 20 Ma due to combined slab and upper plate extensional forcing. At t * = 1, the weak heterogeneity morphology forms a long linear east-west track, extending from the trench to 25 cm (~1625 km) in the east. This feature is extremely thin in the north-south extent (~4 cm, ~260 km). Comparatively, lower north-south migration and less deformation occur when high viscosity heterogeneity is present in the wedge (Fig. 16d -f). Trench-parallel translation rates up to ~3 cm min -1 move residuum south by ~7 cm (~450 km) over 20 Ma. The strong residuum remains robust throughout and trench-normal extension is restrained due to the east-west stress guide.

Discussion
The experiment results produce general aspects of residuum interaction in 4D circulation and specific results relate to the Cascades system. Typically, plume models attribute the Columbia River/Steens flood basalts to processes related to plume head interaction with the lithosphere [e.g., Richards et al., 1989;Takahahshi et al., 1998;Camp and Hanan, 2008;Wolff et al., 2008;Obrebski et al., 2010]. Subsequent eastward trending SRP volcanism is the surface expression for the long-lived conduit tail [e.g., Pierce and Morgan, 2009;Kincaid et al., 2013], and the westward trending HLP volcanic track is due to the deformations and entrainment of remaining plume head modified upper mantle [e.g., Jordan et al., 2004;Kincaid et al., 2013]. Non-plume models invoke plate-driven processes, rather than a mantle plume, to control mantle dynamics and volcanism in the Pacific Northwest.
Subduction-related processes such as slab steepening, rollback, and extension in the overriding plate lead to a pulse of rapid mantle upwelling responsible for the CSFB event [e.g., Faccenna et al., 2010;Liu and Stegman, 2012;Long et al., 2012].  Our models that include back-arc extension in the overriding plate are the most successful in linking experiment results with geological observations of the PNW (e.g., Fig. 17). Interaction between wedge upwelling induced by the extension center and a deforming residuum results in time-progressive melting patterns. In the previous section we resolved that westward age-progressive melting west of the BAE axis, i.e., beneath the HLP, is readily explained in cases with strong or weak residuum present.
However, there is one preferred model that appears capable of upwelling and possible eastward progressing lavas under SRP. Cases with strong heterogeneity located under an overriding plate with BAE result in the largest upwelling east of the extension location, beneath the SRP region (Fig. 14c). High viscosity residuum that resists deformation tends to alter the natural pattern in upwelling along the slab centerline of this system. The resistance to deformation manifests itself as diverted and focused upwelling along the eastern and western edges of the residuum as depicted in schematic models, i.e., region beneath HLP and SRP [e.g., Humphreys et al., 2000;Smith et al., 2009]. This pattern tends to only occur where residuum is present and so the progressive accumulation of this material towards the centerline also focuses the trend of east end and west end upwelling along tracks at the centerline. This is consistent with north-south offset of the CSFB and the HLP/SRP and the centered, There is an interesting dissimilarity between our models and the model proposed in Humphreys et al. [2000]. The mantle residuum in Humphreys et al. [2000] grows in volume over time. Decompression melting occurring along the edges of the viscous heterogeneity creates more refractory residuum along the residuum interface. Through time the feature expands trench-normally (to the east and west) and bimodal melting progresses further along the HLP and SRP. Our models also show upwelling related decompression melting along the eastern and western sides of the high viscosity heterogeneity beneath the HLP and SRP tracks. However, in our simplified fluid experiments, which utilize a thermally altered path of glucose syrup as a mantle residuum analog, the heterogeneity is free to deform but is limited to the same volume throughout the course of the experiment. We can theorize that creating new residuum along the edges in our model would also deflect the upwelling further eastward and westward. This additional model complexity would foreseeably create eastwardprogressing melting beneath the SRP track.

Conclusions
Plate-driven circulation in the subduction system is 3D and time evolving.  Columns 2-3 represent the average plate forcing conditions for each experiment with U D and U T representing the downdip (or longitudinal) sinking and translational rollback speeds of the slab. In the lab, 1 cm min -1 is approximately equivalent to 0.8 cm yr -1 in the mantle. Dip angle ( θ ) is given in degrees from horizontal. Columns 4-6 list the wedge parameters, which include the mantle surface conditions for each experiment (OP: with overriding plate, BAE: with back-arc extension, and blank: free slip mantle surface conditions), mantle wedge fluid temperature (T M ) and corresponding viscosity (μ M ). Columns 7-8 list the initial central residuum location (x i ; y i ) where x i is the distance from the trench and y i is the distance from the slab centerline. For the majority of experiments the residuum is initially located in the center of the wedge (x i = 14; y i = 0), except for the Cascades cases where the residuum location represents the northern wedge location of the CSFB (x i = 13; y i = 6). Residuum material temperature (T R ) and corresponding viscosity (μ R ) are shown in Columns 9-10. Column 11 lists the residuum-wedge viscosity ratio μ R / μ M (μ * ) for each experiment (μ * < 1: weak mantle heterogeneity present in the mantle wedge, μ * > 1: strong mantle heterogeneity present, and μ * = 1: reference cases, no heterogeneity present).

OP$
Le?"behind" as"slab" migrates" westward"   Figure 13. Normalized flow direction (Φ*) and mean instantaneous wedge velocity (U M /U D ) for BAE cases Exp. 27 (weak residuum) and Exp. 18 (strong residuum). The majority of the residuum mass is located in the slab centerline region and velocities in this region exhibit the greatest variation between weak and strong residuum cases.