Abstract

The migration pathways of hydrous fluids in the mantle wedge are influenced by the compaction of the porous mantle matrix, which depends on the matrix permeability, fluid viscosity, and fluid density. Experimental studies show that when fluids are interconnected, the permeability depends on mineral grain size and porosity, the latter of which depends on the amount of fluids introduced into the system (fluid influx). Here, we investigate the role of fluid influx, fluid viscosity, and fluid density in controlling fluid migration in the mantle wedge, using a 2-D numerical model accounting for the effects of grain-size variation and matrix compaction. Our models predict that fluid influx and fluid viscosity are key controls on fluid pathways, while fluid density plays a secondary role. Temperature dependence of fluid viscosity promotes downdip drag of fluids at the base of the forearc mantle toward the subarc region. High fluid influx at postarc depths promotes updip flow near the base of the mantle wedge, guiding the fluids arcward. The model that is applied to northern Cascadia predicts upward fluid migration focused beneath the arc but cannot explain high electrical conductivity observed slightly west of the upward fluid migration. We estimate the amount of hydrous melt that can be produced in the mantle wedge using calculated fluid distributions. Up to a few percent partial melting is predicted in a relatively small region in the core part of the subarc mantle wedge in most subduction settings, including northern Cascadia, and beneath the backarc in old-slab subduction zones.

1. INTRODUCTION

Generation of arc magmas in subduction zones occurs primarily because the addition of water from the dehydrating slab to the overlying mantle lowers the solidus of the mantle (Tatsumi et al., 1986; Grove et al., 2006). Experimental studies and numerical thermo-petrological models predict that the depths of slab dehydration depend on the thermal state of the slab and can vary widely among different subduction zones (Schmidt and Poli, 1998; Hacker et al., 2003; van Keken et al., 2011), potentially promoting magma generation over a wide region. However, present-day arcs are located in a relatively narrow region that lies ∼100 km above the top of the slab regardless of the thermal structure of the subduction zone (England et al., 2004; Syracuse and Abers, 2006; England and Katz, 2010), and it is unclear what controls the location of the arc. What are the pathways of hydrous fluids in the mantle wedge? How do fluids become focused beneath the arc?

Theoretical and experimental studies predict that fluids flow through interconnected pores at grain edges at the pressure and temperature conditions relevant to the ductile part of the mantle wedge (Bulau et al., 1979; von Bargen and Waff, 1986; Mibe et al., 1999), and the governing equations that describe the percolation of fluids through a viscously deformable porous matrix (Sleep 1974; McKenzie 1984; Scott and Stevenson, 1984; Fowler, 1985) have been widely used to model the migration of slab-derived fluids and mantle melts in subduction zones. In this physical model, fluid migration is controlled by solid flow and fluid flow through the solid matrix, the latter of which is driven by the buoyancy of the fluid and the pressure gradients that are induced by compaction and dilation of the viscous matrix over a characteristic length, called the compaction length (Spiegelman, 1993). Fluid flow also depends on fluid mobility, defined here as the ratio of the matrix permeability to fluid viscosity. The former is related to the geometry of the pores network and is proportional to the square of mineral grain size and to a power (generally between 2 and 3) of the fluid fraction in mantle rocks (von Bargen and Waff, 1986; Wark et al., 2003; Simpson et al., 2010; Miller et al., 2014). Based on this model of fluid transport, several hypotheses regarding how fluids are focused beneath the arc have been introduced.

In an earlier study, Spiegelman and McKenzie (1987) proposed the role of dynamic pressure as the main driving force to redirect fluids toward the subarc region. In their model, large dynamic pressure gradients induced by mantle corner flow drove fluids toward the tip of the corner flow. However, the magnitude of pressure gradients was overestimated by the use of a uniform mantle viscosity of 1021 Pa s, which largely exceeds current estimates of temperature-dependent and strain-rate–dependent viscosity thought to better represent mantle rheology (e.g., Kelemen et al., 2003; van Keken et al., 2008).

Wilson et al. (2014) studied the role of matrix compaction in controlling fluid pathways in subductions zones, using a temperature-dependent solid rheology in their fluid migration models. Based on their results, they found two possible mechanisms that focus fluids that are liberated from the subducting slab at postarc depths toward the subarc region, emphasizing the role of compaction in controlling fluid pathways. The first mechanism occurs in the uppermost part of the subducting slab where fluids are produced by dehydration and form a layer of high permeability parallel to the dip of the slab. Resulting compaction pressure gradients further drive nearby fluids into this region, resulting in updip channelized fluid flow through the uppermost part of the slab. Thereby, some of the fluids produced at postarc depth migrate updip and are subsequently introduced into the mantle wedge at subarc depths. The second mechanism occurs in the mantle wedge. The fluids that do not migrate updip are directly released into the mantle wedge. During their ascent through the mantle wedge, their pathways are deflected toward the trench. The deflection is promoted by an upward increase in solid shear viscosity in the shallowest part of the mantle wedge, leading to the formation of decompacting channels that are subparallel to the isotherms. The results showed that the fluid pathways depend strongly on the compaction length, over which fluids feel the effects of matrix compaction, and that the two above mechanisms are only effective when compaction lengths are sufficiently high. Because compaction depends on fluid mobility, the authors suggested that spatial variation of grain size is another possible factor that contributes to focusing of fluids beneath the arc because it affects the mantle permeability.

Using simplified models of porous flow, in which the compaction of the matrix was neglected, several studies address the migration of hydrous fluids in the mantle wedge (Iwamori, 1998, 2000; Arcay et al., 2005; Cagnioncle et al., 2007; Hebert and Montesi, 2013; Angiboust et al., 2012; Wada and Behn 2015). Cagnioncle et al. (2007) illustrated the influence of grain size on fluid flow in the mantle wedge. Assuming a uniform grain size, the authors showed that relatively small grain size could lead to trapping of fluids at the base of the mantle wedge and their downdip transport by mantle outflow. Later, Wada et al. (2011) estimated the grain-size distribution in the mantle wedge using the “paleowatmetter” model of Austin and Evans (2007). Grain size was shown to increase from a few µm near the tip of the corner flow to a few cm in the hottest part of the mantle wedge. The calculated grain-size distribution was adapted by Wada and Behn (2015) in a simplified porous flow model to test its influence on fluid migration. The results showed that small grain size along the base of the mantle wedge at shallow depth causes entrapment of fluids in the downgoing mantle, providing a mechanism for downdip transport of fluids from forearc depths to the subarc region. In their model, however, the effect of matrix compaction was not incorporated.

In our previous study, we developed numerical models of fluid migration in the mantle wedge accounting for both grain-size distribution and mantle compaction (Cerpa et al., 2017). The modeling results showed that spatial variations of grain size and solid shear viscosity both serve to focus fluids beneath the arc, and that their contributions as fluid-focusing mechanisms vary spatially in the mantle wedge. Grain-size variation plays a key role at the base of the mantle wedge near the tip of the corner flow and focuses fluids released beneath the forearc toward the subarc region for all modeled subduction zones as long as fluids are introduced beneath the forearc. The pathways of fluid released at postarc depths in cold-slab subduction zones are deflected toward the arc at shallow depths due to gradients of compaction pressure and advection by mantle inflow. In the study, the pattern of fluid release from the slab was calculated based on the thermal state of the slab, but only a fraction of the total fluid released from the slab was introduced into the mantle wedge such that upward fluid migration occurs in the model. Fluid influx is an important parameter that determines the fluid fraction (i.e., porosity) at the base of the mantle wedge and thus partly controls the permeability in this region and subsequent evolution of fluid migration pathways in the overlying mantle. Fluid viscosity and fluid density vary by several orders of magnitude with temperature and composition (mainly the water content) at relatively high pressure (e.g., Audétat and Keppler, 2004) and can also affect fluid migration, yet their effects were not quantified in the previous study.

The present work focuses on quantifying the effects of fluid influx and fluid properties on the pattern of fluid migration through the application of the fluid migration model developed by Cerpa et al. (2017). In the following, we first briefly summarize the results of theoretical and laboratory studies, from which our model parameters are derived, and the modeling approach. We then present a series of modeling results that illustrate the effects of fluid influx, fluid viscosity, and fluid density. In the last three simulations presented in this study, the subduction parameters are chosen to represent the Northern Cascades subduction zone, where the distribution of fluids in the mantle wedge has been inferred from magnetotelluric (MT) observations. Finally, we discuss our results with an emphasis on the implications for hydrous melting in the mantle wedge.

2. PERMEABILITY AND FLUID PROPERTIES IN THE MANTLE WEDGE

2.1. Permeability of Mantle Rocks

Fluid migration by porous flow occurs if the fluids form an interconnected network of channels or pores within the hosting solid rock. Indeed, an intergranular fluid phase can be distributed along grain edges forming the channels or at the grain corners forming nodes of fluid, depending on which configuration minimizes the free energy at grain boundaries (Bulau et al., 1979). The dihedral angle θ, that is, the angle between two fluid-solid interfaces at a corner of a fluid-filled pore, controls the topology and interconnectivity of the fluid phase in texturally equilibrated rocks. When θ < 60°, channels are interconnected, and when θ > 60°, the channels are closed off, and fluid is trapped in isolated pockets at grain corners, preventing fluid migration (Bulau et al., 1979; Watson, 1982). However, von Bargen and Waff (1986) computed dihedral angles for single-phase isotropic tetrakaidekahedral crystals of constant size and found that the transition from open to closed channels (“pinch-off” criterion) only applies at relatively low fluid fractions (<<0.01), and at relatively high porosities (>>0.01), grain boundaries may become wetted (Waff and Faul, 1992; Miller et al., 2014), particularly when the bulk undergoes high shear stress conditions (Hier-Majumder and Kohlstedt, 2006). At the base of the mantle wedge, a relatively high fluid fraction and a high shear stress condition are expected, likely promoting interconnectivity.

Laboratory measurements of the dihedral angle in synthetic rock samples hosting CO2-H2O fluids at 1000 °C and 1 GPa were carried out by Watson and Brenan (1987). Their results showed high dihedral angles (above 70°), suggesting that hydrous fluids could not form an interconnected network in the mantle at relatively low porosities. However, in a subsequent study, Mibe et al. (1998) measured the dihedral angle in aqueous fluid-forsterite systems at 1000 °C and 3–5 GPa and reported a mean dihedral angle of 40°. They emphasized that θ may drop below 60° at 2 GPa as a result of changes in pressure-dependent olivine solubility that affects the composition of hydrous fluids. Therefore, at the conditions of the flowing part of the mantle wedge, fluids likely occupy interconnected pores and migrate by porous flow (Mibe et al., 1999).

At textural equilibrium, the relationship between permeability, fluid fraction, and grain-size follows a power law of the form of 
graphic
where b is the grain size diameter or channel spacing, ϕ is fluid volume fraction or porosity, n is a porosity exponent, and C is a geometrical factor (McKenzie, 1984; von Bargen and Waff, 1986; Turcotte and Schubert, 2014). In this study, we use C = 270 (Wark et al., 2003).

In the creeping regime, grain size can evolve through grain growth and dynamic recrystallization (Karato, 1984). Grain growth occurs by grain boundary migration, minimizing the interfacial free energy between grains (Atkinson, 1988; Evans et al., 2001). Under applied stresses, dynamic recrystallization occurs by reorganizations of dislocations to form new subgrain boundaries and/or by migration of existing grain boundaries to form new grains.

Austin and Evans (2007) derived a model for the evolution of grain size called the “paleowattmetter” model, in which the rate of mean grain-size reduction through dynamic recrystallization is related to the mechanical work done to the grains. In the model, dynamic recrystallization and grain growth are assumed independent of each other, and when their rates are balanced, the mean grain size reaches its steady state (Behn et al., 2009). Wada et al. (2011) coupled the “paleowattmetter” model to thermomechanical models of subduction zones to estimate the steady-state grain-size distribution in the mantle wedge. They found that the grain size increases from a few µm near the tip of the corner flow to a few cm in the hottest part of the flowing mantle wedge and that this downdip increase in grain size by approximately two orders of magnitude occurs regardless of the subduction parameters. The model likely overestimates the rate of downdip increase in grain size with depth at the base of the mantle wedge where the effect of advection, which is neglected in the grain-size calculation, becomes relatively important (Wada et al., 2011), and therefore the fine grain size may persist to even greater depths than the model predicts. Furthermore, the “paleowattmeter” model is based on experiments performed with single-phase (olivine) systems and neglects the effects of a secondary phase on grain growth. The presence of a secondary phase, such as pyroxene, likely limits the growth of olivine grains to a few cm (Evans et al., 2001, and references therein). Nonetheless, Wada et al. (2011) tested the effects of imposing an upper limit of their calculated grain size to a few cm and found little effect of this value on the grain-size distribution in the regions where the grain size is smaller than a few cm.

The effects of grain-size distribution on fluid migration have been studied by Wada and Behn (2015) and Cerpa et al. (2017), who showed that fluids become trapped in small grain size at shallow depths and become dragged downdip by solid flow until grain size becomes large enough to allow upward fluids migration, potentially redirecting shallow fluids toward the subarc region. While these studies focused on the effect of grain size on permeability, the role of fluid influx at the base of the mantle wedge was not investigated. The imposed fluid influx controls the fluid fraction and thus permeability at the base of the mantle wedge. Thus, in this study, we focus on the role of fluid fraction at the base of the mantle wedge in controlling the evolution of fluid migration.

The parameter n from Equation (1) that relates permeability to fluid fraction has been estimated experimentally and numerically for partially molten rocks. There has been little experimental work done on the quantification of n for a hydrous-fluid–bearing rock. However, the theoretical models used for the prediction of n are based primarily on geometrical considerations, and therefore the results have been applied to models of water-rich fluids migration. A similar reasoning applies to the geometrical factor C.

Estimations of n range from ∼2 for idealized isotropic systems of uniform grain size (Frank, 1968; Maaløe and Scheie, 1982; von Bargen and Waff, 1986; Simpson et al., 2010) to ∼3 for more complex systems (anisotropic, non-uniform grain sizes) (Wark et al., 2003; Connolly et al., 2009). More recently, Miller et al., (2014) conducted laboratory estimations of permeability in synthetic partially molten rocks using X-ray microtomography and have inferred an exponent of n = 2.6 from their data set. These authors observed two topological regimes that could be represented by two end members: n = 2, when melt resides principally along channels at relatively low porosities and n = 3, when melt tends to form films at grain boundaries at relatively high porosities (see also Rudge, 2018). As discussed by the authors, their inferred value of n may represent the coexistence of both regimes in different sub-volumes of their samples. In our study, we are primarily interested in melt migration from the bottom of the mantle to the region of flux melting where a relatively high fluid fraction is expected. Furthermore, the mantle wedge corner flow is a high-deformation environment, which is likely to promote wetting of grain boundaries (Hier-Majumder and Kohlstedt, 2006). Hence, in our fluid-flow models, we assume n = 3.

The fluid influx at the base of the mantle wedge controls the overall porosity in the system. To assess the effects of fluid fraction on permeability and on fluid-flow pathways, we systematically vary the fraction of the thermopetrologically predicted slab-derived fluids that are allowed to flux into the mantle wedge in our model.

2.2. Fluid Viscosity and Density

A series of dehydration reactions occurs in the downgoing slab, producing a hydrous fluid that can migrate into the mantle wedge. The hydrous fluid likely dissolves silicate minerals as it percolates the subducted lithosphere. The solubility of silicates in hydrous fluids increases with pressure and temperature (Manning, 2004, and references therein), and thus the distribution of dissolved silicates in the fluids that migrate into the mantle wedge depends on the thermal state of the slab and history of the fluid. Fluids that migrate into the hot part of the mantle wedge trigger melting of the mantle (e.g., Grove et al., 2006). The resulting silicate-rich melts contain some amount of water as inferred from the H2O contents found in arc lavas (Wallace, 2005, and references therein), but the exact H2O contents in arc magmas are still debated (Plank et al., 2013). Hence, from the slab to the overlying crust, the composition of the migrating fluids can vary significantly, resulting in spatial variations in the physical properties of the fluids. Of those, viscosity and density affect the mobility and the buoyancy of fluids, respectively, and are critical parameters to fluid migration.

Fluid viscosity in the mantle wedge is expected to vary by several orders of magnitude, spanning from values similar to pure water 10−4 Pa s near the base of the mantle wedge to those of highly silicic melts 1011 Pa s as they reach the overlying crust (Audétat and Keppler, 2004; Giordano et al., 2008; Hack and Thompson, 2011). Audétat and Keppler (2004) measured the viscosities of H2O-rich fluids (>20 wt%) with varying water and dissolved silica contents at high temperatures (600–950 °C) and high pressures (1–2 GPa). They found that fluid viscosity is insensitive to pressure variation but is strongly dependent on temperature and composition and proposed an empirical relationship for the viscosity (µ) of fluids as a function of dissolved silica concentration (csil) (in wt%) and temperature (T) (in K) in subduction zones: 
graphic

The authors suggest that their model reproduces the viscosities of fluids that contain up 80 wt% of dissolved silicate in experiments. Their model predicts that fluid viscosity varies from 10−4 Pa s (e.g., viscosity of pure water) for csil = 0 wt% to 10 Pa s for csil = 80 wt% at T = 700°C (Fig. 1). In this study, we test the role of fluid viscosity in controlling fluid migration pathways by applying a range of constant fluid viscosities in the range of 0.001–10 Pa s in our fluid migration models, which corresponds approximately to a temperature range of 700–1300 °C and dissolved silica concentrations of >20 wt%.

In contrast to fluid viscosity, fluid density is not expected to vary by more than a factor of two in the mantle wedge. We will illustrate the effects of density through three models with different density values.

3. MODELING APPROACH

We use the modeling approach of Cerpa et al. (2017), which incorporates the effects of grain-size distribution and mantle compaction on fluid migration. Given the relatively low porosities that are calculated in our models, it is reasonable to assume that the effects of the fluid phase on the solid-state mantle flow are negligibly small (Wilson et al., 2014). Based on this assumption, we neglect the effect of compaction on the divergence of solid velocity field (e.g., Dymkova and Gerya, 2013), the effects of porosity on the buoyancy of the mantle, and the weakening of the solid due to the presence of a free-fluid phase (Hirth and Kohlstedt, 2003. The advection of heat by the fluid phase (Rees-Jones et al., 2018) is also neglected in the temperature calculations.

The governing equations are the same as in Cerpa et al. (2017) and are briefly summarized in this section. Parameters with a prime are non-dimensional, while other parameters are dimensional. The subscript 0 indicates that the parameter is a reference parameter used for the non-dimensionalization of the equations.

3.1. Governing Equations

3.1.1. Solid Flow

We assume that the solid-state mantle flow has reached a steady state. The governing equations for the solid flow are: 
graphic
 
graphic
 
graphic
where forumla is mantle (solid) velocity vector, p is dynamic pressure, T is temperature, κ is thermal diffusivity, Q is radiogenic heat source, and ⋵ is strain rate tensor defined as 
graphic

The mantle (solid) shear viscosity η accounts for both diffusion creep and dislocation creep (see Appendix A). The above governing equations are exactly the same as those in Cerpa et al. (2017) except that we include radiogenic heat in the heat transfer equation.

As in Cerpa et al. (2017), the thermomechanical model is coupled to the paleowattmeter model (Austin and Evans, 2007) to calculate the grain-size distribution in the mantle wedge. Our model assumes that the timescale of grain-size evolution is small compared to the timescale of the fluid migration, and thus we neglect advection in the calculation of the steady-state grain size (Behn et al., 2009; Wada et al., 2011). The parameters used in the calculations are given in Appendix A.

3.1.2. Fluid Flow

In the fluid-flow calculations, we neglect mass exchanges between phases, and we do not account for the effects of the dynamic pressure due to mantle shear. The non-dimensional governing equations for the fluid flow are thus 
graphic
 
graphic
where ϕ is porosity or fluid fraction, Ƥ is compaction pressure, ζ and K are bulk viscosity and permeability of the solid, respectively, µ is fluid viscosity, ĝ is the unit vector in the direction of gravity, forumla are fluid and solid reference velocities, respectively. h0 is a reference length scale and δ0 is a reference compaction length. The relationships among the reference parameters that are important in our study are given in Appendix B along with their values (see also Cerpa et al., 2017). In particular, in our models the reference compaction length is defined as 
graphic
where Δρg is fluid buoyancy, η0, K0 and µ0 are reference solid viscosity, permeability, and fluid viscosity, respectively. Ψ0 is a reference flux used for the non-dimensionalization of the fluid influx boundary condition as described below. As described by Wilson et al. (2014), compaction lengths govern the mode of fluid flow. Low compaction lengths favor a fluid transport by solid advection, while high compaction lengths promote a buoyancy-driven flow and enhance the effects of compaction pressure. In our models, the compaction length varies spatially and temporally, and the non-dimensional compaction length is defined as 
graphic
where permeability (K′) and bulk solid viscosity (ζ′) are 
graphic
 
graphic

We regularize the above two parameters to avoid numerical singularities (Wilson et al., 2014; Cerpa et al., 2017) and express regularized values as K′ and ζ′.

The fluid velocity, υ′f , is the vectorial sum of the solid velocity and the fluid velocity relative to the mantle matrix based on Darcy’s law: 
graphic

3.2. Model Setup

The model setup for solid-state mantle flow is similar to that used in Cerpa et al. (2017), except that the Cascadia model accounts for internal radiogenic heat in the crust (Wada and Wang, 2009).

In the solid-state mantle flow model, the subduction velocity (υsub) is prescribed at the top and the bottom of the subducting slab. We assume the slab and the overlying mantle are decoupled at shallow depths and become fully coupled at a 75 km depth (hereafter referred to as the maximum decoupling depth [MDD]) to be consistent with previous studies (e.g., Wada and Wang, 2009). The rigid overriding crust is fixed. The geotherm applied on the trench-side vertical boundary of the model is calculated using the GDH1 plate cooling model, which assumes a 95-km-thick slab (Stein and Stein, 1992). The geotherm imposed on the shallow part of the backarc-side boundary is calculated using a one-dimensional conductive heat model with an average surface heat flux of 80 µW/m−2. As stated above, internal heat production is included in Cascadia models (Table 1). On the deeper part of the backarc side boundary, down to the approximate depth of the transition between mantle inflow and mantle outflow, we impose a geotherm calculated by a potential temperature of 1350 °C and an adiabatic gradient of 0.3 °C km−1. The surface temperature (Tsur) is 0 °C, and slab bottom temperature (Tbot) is 1450 °C.

Following Wada et al. (2012) and Cerpa et al. (2017), we estimate the pattern of fluid release from the slab by first computing the temperature field and then computing the water produced by dehydration reactions, using tabulated phase diagrams from the thermodynamic calculation code Perple_X (Connolly, 2009). We assume a hydrated layer at the uppermost part of the subducting plate arriving at the trench consisting of a 0.5-km-thick volcanic layer, a 1.5-km-thick layer of dikes, a 5-km-thick gabbros layer, and a 4-km-thick peridotitic uppermost mantle layer with initial water contents of 2.6, 1.8, 0.8, and 2.0 wt%, respectively. The uppermost part of the slab is divided into 1-km-wide columns, each consisting of 500-m-tall cells. Assuming vertical fluid flow within the slab, the integration over each column of the water release per cell provides an estimation of the volumetric flux QH2O of fluids released at the top of the slab as a function of depth (z). Because our study focuses on fluid migration in the creeping mantle wedge, we exclude the “cold nose,” that is, the cold corner tip between the subducting slab and the overlying crust that does not deform, from our calculations, and we estimate fluid release from the slab only below the MDD. In the calculation, the effects of updip migration (Wilson et al., 2014) and rehydration reactions (Wada et al., 2012) within the slab are neglected. The fluid influx to be applied in the model (Q (z)) is a product of the calculated QH2O (z) and a tuning parameter ω, which controls the fraction of slab-derived fluids that are introduced into the mantle wedge to account for the fact that all fluids released from the slab do not migrate instantaneously into the mantle wedge. In Cerpa et al. (2017), we applied a few different fluid influx patterns QH2O (z) to test their effects on fluid pathways in the mantle wedge. The tuning parameter value that controls the magnitude of fluid influx was fixed for all the models and was chosen as the minimum value that allowed upward migration in all the models. Here, we investigate the influence of the tuning parameter ω.

In the fluid-flow calculations, Equations (7) and (8) are not solved in the slab nor in the overriding crust. The non-dimensional volumetric flux Q′(z) = Q(z)/Ψ0 is imposed as a boundary condition at the base of the mantle wedge. The top of the fluid domain through which fluids are allowed to escape is set at a 45 km depth. Above this depth, fluid migration mechanisms involving mechanical failure of rocks may play a role (Havlin et al., 2013; Keller et al., 2013) but are neglected in our models. As described in Cerpa et al. (2017), in the fluid calculations, we use the mantle shear viscosity from the solid calculation to which we impose a cap value of 1021 Pa s. Similarly, we impose a lower bound of 10 µm and an upper bound of 2 cm to the grain size in the fluid-flow models.

All calculations are performed using the software package TerraFERMA (Transparent Finite Element Rapid Model Assembler) (Wilson et al., 2017).

4. EFFECTS OF FLUID INFLUX AND FLUID PROPERTIES ON FLUID MIGRATION

We first present models with a simple subduction zone geometry with a constant slab dip of 30° and a 30-km-thick overriding crust. For the calculation of the thermal and solid-state mantle flow fields, we impose a subduction rate of 4 cm/yr and a slab age of 50 Ma. In the model (S1), full coupling between the slab and the overlying mantle downdip of the MDD induces a corner flow that replenishes the mantle wedge with the hot mantle from the backarc region (Fig. 2A; data used to create figures in Supplemental Material1). The model-predicted distributions of mantle wedge flow, shear viscosity, and grain size have been described in Cerpa et al. (2017). The key observations (Figs. 2B and 2C) that are important to fluid migration include (1) low shear viscosity (∼1018 Pa s) along the line separating the shallow incoming mantle flow toward the trench from the deeper downgoing mantle flow due to relatively high shear stresses and temperature; (2) moderate viscosity (1019–1020 Pa s) in the inflow region where the strain rate is relatively low; (3) increasing viscosity (>1020 Pa s) toward the surface above the mantle inflow region because of an upward decrease in temperature; and (4) downdip and arcward increase in grain size from a few µm near the tip of the corner flow to a few cm in the hottest part of the mantle wedge.

The thermal structure of model S1 leads to release of fluids from the slab over two distinct depth ranges (Fig. 3). The first peak in the fluid release between 75 and 100 km depths is associated with the dehydration of the upper crust while the second peak between 135 and 150 km depths is associated with the dehydration of the subducting mantle.

We use the above solid shear viscosity structure, grain-size distribution, and pattern of fluid release in the fluid-flow calculations with a simple slab geometry (models F1–11) as described in Section 3. In the calculations, we apply uniform fluid properties unless otherwise stated.

4.1. Effects of Fluid Influx

We use a model with a fluid viscosity of 0.1 Pa s and a tuning parameter of 0.5% as a reference model (model F1). The model is identical to model 50E in Cerpa et al. (2017), and therefore we only summarize the main features. Fluids introduced during the first peak dehydration are initially dragged downdip (Fig. 4A) because of the small grain size near the MDD and the presence of a layer of relatively high solid shear viscosity that develops just above the slab. As grain size increases downdip, the mantle permeability and the compaction length increase, and fluids start to migrate upwards at subarc depths (Figs. 4A and 4B). Fluids released by the second peak dehydration migrate largely upward because of the relatively high grain size in the deep mantle wedge (Figs. 4A–4C). In the mantle inflow region, fluid pathways are deflected trenchward under the effects of the compaction pressure gradients and the solid advection (Fig. 4C). In the shallowest part of the mantle wedge where solid flow vanishes, fluids tend to pond due to the vertical increase in solid shear viscosity and decrease in grain size. In our fluid-flow models, fluids initially migrate through a dry mantle. As we are rather interested in mature fluid migration pathways (Figs. 4C and 4D), we use the time-averaged porosity field (Fig. 4E) to represent the “long-term” behavior of the fluids.

We now consider a model with a higher tuning parameter of 5% (model F2), i.e., a larger flux of fluids introduced at the base of the mantle wedge. All other parameters remain identical to model F1. In this model, fluids from the first peak dehydration are first dragged downdip as in model F1 and then start to migrate upward at subarc depths (Fig. 5A). However, the distance over which these fluids are dragged downdip is shorter compared to that in model F1. Some fluids from the second peak dehydration migrate immediately upward as in model F1, but there are also some fluids that first migrate updip along the base of the mantle wedge (Fig. 5A). Above the second peak dehydration, solid viscosity is relatively high, and compaction pressure gradients guide fluids to migrate updip through high permeability layers oriented parallel to the dip of the slab. Those fluids start to move upwards at a ∼130 km depth, but during their ascent, some of them become entrained and dragged downdip by the outflowing mantle, merging back with the rising fluids from the second peak dehydration (Fig. 5B). This complex pattern equilibrates progressively, eventually leading to the formation of two main vertical fluid pathways separated by a horizontal distance of ∼30 km from each other (Figs. 5C and 5D). In model F2, upward migration of fluids is less affected by mantle inflow and variations in solid shear viscosity in flowing mantle wedge than in model F1 owing to higher permeability due to increased fluid influx. In the cold shallowest mantle wedge, fluids tend to pond as in model F1. The time-averaged porosity field (Fig. 5E) shows wider fluid pathways than in model F1. In the following, only time-averaged porosity is discussed.

We test the effects of a high tuning parameter (50%) in model F3 (Fig. 6), while all other parameters remain identical to models F1 and F2. In this model, the distance over which fluids from the first peak dehydration are dragged is even shorter compared to that in models F1 and F2. Conversely, the distance of updip migration of some fluids from the second peak dehydration is longer than that in model F2, resulting in two vertical fluid pathways about a 50 km distance apart. In model F3, solid advection plays a modest role in controlling fluid migration.

4.2. Effects of Fluid Viscosity

4.2.1. Uniform Fluid Viscosity

Relative to the reference model F1, we decrease and increase fluid viscosity by two orders of magnitude to 0.001 in models F4 and 10 Pa s in F5, respectively (Fig. 6). All other parameters are identical to model F1. The fluid viscosities in models F4, F1, and F5 correspond to dissolved silicates concentrations (csil) of 20.5, 54.6, and 88.8 wt%, respectively, at a temperature of 800 °C based on the empirical law for aqueous fluid viscosity of Audétat and Keppler (2004) (Fig. 1).

Decreasing fluid viscosity promotes greater upward migration from the first peak dehydration and updip migration of fluids from the second peak dehydration. Conversely, increasing fluid viscosity by two orders of magnitude leads to the entrapment of fluids at the base of the mantle wedge and their downdip drag by the mantle flow.

Although the absolute values of fluid fraction in model F4 differs from that in model F2, the two models show identical fluid migration pathways, indicating that decreasing fluid viscosity by two orders of magnitude has the same effect as increasing fluid influx by one order of magnitude. In other words, models with the same ω2/µ produce the same fluid pathways although the distributions of fluid fraction are different due to the difference in the amount of fluids introduced into the system. This scaling relationship in controlling compaction-driven flow is consistent with the relationship described in Equation (9) for the reference compaction length, fluid flux, and fluid viscosity; the compaction length is scaled with Ψ020, where Ψ0 plays an equivalent role to ω as the latter controls fluid influx.

We include two additional models F6 and F7 to demonstrate this scaling relationship. In these models, ω and/or µ take different values from those in models F1–F5, but ω2/µ in F6 is identical to that of model F1, and that in F7 is identical to those in model F2 and F4, and models that share the same value show identical fluid migration paths. The tuning parameter ω controls fluid influx and thus the porosity and permeability at the base of the mantle wedge, affecting fluid mobility and compaction length. As a consequence, for a given grain-size distribution, the impact of matrix compaction is mainly controlled by ω2/µ.

Models with ω2/µ < 2.5 × 10−6 Pa−1 s−1 (model F5) lead to entrapment of fluids at the base of the mantle wedge and little upward fluid migration. In models with ω2/µ < 2.5 × 10−4 Pa−1 s−1, downdip drag is still prominent for fluids from the first peak dehydration, but upward fluid migration develops. Models with ω2/µ < 2.5 × 10−2 Pa−1 s−1 show a lesser effect of solid advection and eventually promote some updip fluid migration of fluids from the second peak dehydration at the base of the mantle wedge.

4.2.2. Temperature-Dependent Fluid Viscosity

In nature, fluid viscosity is expected to vary spatially because of its dependence on fluid composition, pressure, and temperature. Variation in the composition of hydrous fluids that are traveling through the mantle wedge depends largely on silicate solubility, which increases with both the temperature and the pressure (e.g., Newton and Manning, 2002). As fluids travel to shallower depths from the base of the mantle wedge to its hot core, pressure decreases but temperature increases. Given the competing effects of temperature and pressure on silica solubility, the change in fluid composition may be relatively small. The direct effect of pressure on fluid viscosity is expected to be also relatively small, whereas temperature has a large impact on fluid viscosity as fluids travel between the relatively cold base and the hot core of the mantle wedge (Audétat and Keppler, 2004). Here, we investigate the effects of temperature-dependent fluid viscosity on fluid migration.

Two end-member models are studied: model F8 with relatively low silicate content csil = 30 wt% and model F9 with relatively high silicate content csil = 70 wt% (Fig. 7). In both models, ω = 5%. Near the “cold nose,” fluid viscosity is expected to be relatively high, while it is expected to be relatively low in the flowing part of the mantle wedge (Figs. 7A–7C). The relatively cold temperature in the shallow mantle wedge beneath the overlying crust and immediately above the 50 Ma slab also results in relatively high fluid viscosity. Fluids are thus expected to be less mobile at the base of the mantle wedge than in its hottest part.

In model F8 with csil = 30 wt%, fluid viscosity at the base of the mantle wedge is 0.01–0.1 Pa s, similar to that in model F2, and it becomes lower (∼10−3 Pa s) in the hot part of the mantle wedge, resulting in a higher fluid mobility in this region compared to model F2. Therefore, model F8 exhibits some downdip drag of fluids from the first peak dehydration at the base of the mantle wedge over a distance similar to that in model F2. The two vertical fluid pathways beneath the backarc are ∼40 km apart and are spaced closer than those in model F2. This indicates that fluids from the second peak dehydration travel farther updip compared to model F2. This is consistent with the higher mobility of fluids expected in the model with temperature-dependent viscosity as they migrate away from the base of the mantle wedge. Fluid pathways at shallow depths are narrower in model F8 than in model F2 because buoyancy-driven upward fluid flow becomes more dominant than fluid migration due to mantle inflow and compaction-driven flow.

In model F9 with csil = 70 wt%, fluid viscosity in the hottest part of the mantle wedge is predicted to be comparable (∼0.1 Pa s) to that of model F2 but higher (10–100 Pa s) at the base of the mantle wedge. As a consequence, fluids from the first peak dehydration become dragged downdip farther than that in model F2, and updip migration above the second peak dehydration does not occur. At the base of the mantle wedge, the effects of solid advection and spatial variations in solid shear viscosity are similar to that in model F2.

Models F8 and F9 show that a temperature-dependent viscosity promotes lateral transport of fluids away from the trench at the base of the mantle wedge. This results in more focusing of fluids that are released beneath the forearc region toward the subarc region. On the contrary, a temperature-dependent viscosity reduces trenchward deflection and subsequent spread of upward fluid migration pathways by mantle inflow at shallow depths. However, as discussed by Cerpa et al. (2017), fluid migration pathways in and above the hottest part of the mantle wedge may differ from our predictions since melting may produce more viscous and less dense fluids.

4.3. Effects of Fluid Density

To test the effects of fluid buoyancy, we reduce the density contrast between the fluid and the solid from the reference value of 2000 kg m−3 to 1000 kg m−3 in model F10 and to 200 kg m−3 in model F11 (Fig. 8). Because Δρ is used for the non-dimensionalization of governing equations, we modify the reference flux Ψ0 and thus the tuning parameter ω accordingly (see Table 2) so that the reference compaction length and boundary conditions in models F10 and F11 are identical to those in model F2 (see Appendix B).

Decreasing the density contrast by a factor of two (model F10) compared to model F2 generates fluid pathways that are similar to those obtained in models F1 and F6. Decreasing density contrast even further, by one order of magnitude (model F11), does have a noticeable effect on fluid pathways, which is comparable to the effect of increasing fluid viscosity by two orders of magnitude (model F5). However, such small density contrast between the fluids and the mantle rock is unlikely, particularly below the core part of the mantle wedge. Therefore, the effect of fluid density variation on fluid migration pathways from the base to the core part of the mantle wedge may be small compared to the effects of variations in fluid viscosity.

5. FLUID PATHWAYS IN THE NORTHERN CASCADIA SUBDUCTION ZONE

The Cascades subduction zone (39°N–53°N) is formed by the subduction of the Juan de Fuca plate beneath the North American plate and is often considered as an end-member warm-slab subduction zone because of the relatively young age (ca. <20 Ma) of the subducting plate (Kirby et al., 1991; Peacock and Wang, 1999; Hacker et al., 2003). Thermal models for Cascadia predict that the subducting Juan de Fuca slab releases most of its water before it reaches the subarc depths (Wada and Wang, 2009; Syracuse et al., 2010; van Keken et al., 2011). However, geochemical and geophysical observations indicate the addition of slab-derived fluids to the subarc mantle, leading to the generation of arc magmas beneath the Cascades.

Based on geochemical analyses of primitive basalts, Schmidt et al. (2008) divided the Cascades margin into four different segments (North, Columbia, Central, and South segments). In particular, the composition of the most primitive calk-alkaline basalts in the Columbia segment between Mount Rainier (46.8°N, 121°W) and Mount Jefferson exhibits high 87Sr/86Sr and Ba/Ce ratios compared to other basalts from the Cascades and may reflect the effects of flux melting (Schmidt et al., 2008, and references therein). Although less pronounced in the Columbia segment compared to the southern part of the Cascades, relatively high ratios of fluid-mobile to fluid-immobile trace elements indicate the addition of slab-derived fluids. There is also some evidence of dry melting in the backarc region ∼46°N (Leeman et al., 2005).

Magnetotellurics (MT) has been used to infer the distribution of fluids in subduction zones, given the sensitivity of electrical conductivity (resistivity) to the presence of fluids and melts (Ni et al., 2011; Pommier and Garnero, 2014). The electrical resistivity model for the region beneath Mount Rainier of the central Cascades shows a high-conductivity anomaly originating just above the slab at ∼80 km depth and extending both upward and downdip with decreasing values in those directions (McGary et al., 2014). Those anomalies were interpreted as indicators of the presence of slab-derived fluids and melts that are generated by flux melting. Another electrical resistivity model for the same region also indicates a deeper high-conductivity anomaly above the slab situated beneath the backarc region at 150–200 km depths (Wannamaker et al., 2014), which is interpreted as a region of decompression melting of upwelling mantle, consistent with geochemical studies in this region (e.g., Leeman et al., 2005). Both MT models of McGary et al. (2014) and Wannamaker et al. (2014) predict the presence of fluids in a region that extends from the base of the mantle wedge to the overlying crust and is deflected toward the trench at shallow depths. This feature was suggested by McGary et al. (2014) to represent the pathways of buoyant diapirs deflected by mantle inflow. The MT results, however, do not provide information on the direction of fluid migration. Quantitative models of fluid flow in the mantle wedge will therefore help to improve the interpretation of MT observations and other geophysical observations, such as seismic velocity and attenuation structures. Here, as a first step, we develop fluid migration models for Cascadia along an E-W profile at 47°N (Fig. 9) to complement the existing MT observations and test the effects of ω2/µ on fluid migration and melt generation.

For the model construction, we adopt the slab geometry model for Cascadia from McCrory et al. (2012). The model extends only to a 95 km depth, and we extend the slab to the bottom of our model, assuming a constant slab dip of 30.3°. In the solid flow calculation, we apply a subduction rate of 3.9 cm/yr, which corresponds to a profile-parallel component of the subduction velocity of the Juan de Fuca plate (4.5 cm/yr) (Wada and Wang, 2009). The subducting plate age at the trench is 7 Ma (Syracuse et al., 2010), and the thickness of the overriding crust is 40 km. Other model parameters for the solid flow calculations are identical to those listed in Table 1.

Compared to the reference model S1 with a 50 Ma slab, the slab in the Cascadia model is much warmer. For instance, the 600 °C isotherm is advected to ∼70 km depth in the Cascadia model (Fig. 10A), while in model S1, the same isotherm extends beyond the bottom of the model domain at 160 km depth. Overall, the grain-size distribution is similar to that in model S1, but grain size is slightly larger particularly at the base of the mantle wedge, owing to the warmer slab thermal state (Fig. 10B). Similarly, the mantle wedge viscosity structure in model SC1 (Fig. 10C) is similar to that in model S1, but the viscosity is slightly lower in model SC1.

The volumetric flux of fluids from the top of the slab in model SC1 is calculated using the same approach described above (Fig. 10D). Because of the high temperatures within the slab, dehydration of the subducting crust occurs beneath the “cold-nose,” but a peak dehydration of the subducting mantle occurs between 75 and 95 km depths, consistent with earlier petrological modeling results for Cascadia (Wada and Wang, 2009; van Keken et al., 2011).

We present three fluid migration models for Cascadia (models FC11–FC13; Fig. 11). In model FC11 and FC13, ω is 5 and 50%, respectively, and µ is 0.1 and 10 Pa s, respectively, yielding ω2/µ of 2.5 × 10−2 Pa−1 s−1. In model F12, ω is 5, and µ is 10 Pa s, yielding ω2/µ of 2.5 × 10−4 Pa−1 s−1. As discussed above, the exact values of ω and µ are not critical to the fluid migration pathways, and the value of ω2/µ depicts the fluid migration paths. In all three models, a constant density contrast of 2000 kg m–3 is assumed. In models FC11 and FC13, fluids immediately flow upward at subarc depths with little initial downdip transport along the base of the mantle wedge. During their ascent, fluids are little affected by solid advection and spatial variations in solid shear viscosity in the inflow region. In model FC12, fluids are transported downdip and then migrate upward at postarc depths (∼90–100 km depths), and some fluids are dragged farther downdip. The deflection of fluids toward the trench by the incoming mantle flow is more pronounced in model FC12. All three models predict some focusing of fluids near the subarc region.

6. DISCUSSION

6.1. Physical Factors Influencing Fluid Pathways in the Mantle Wedge

Wilson et al. (2014) investigated the influence of mantle compaction on fluid migration pathways in subduction zones. In their model, grain size was assumed constant, and mantle permeability varied only with fluid fraction. Thereby, the compaction length was controlled by the reference fluid mobility (i.e., the ratio of matrix permeability to fluid viscosity K00). They showed that, with relatively high fluid mobility (K00 ∼10−7 m3 s kg−1) and thus a high compaction length (δ0 ∼ 20 km), fluid migration paths are strongly influenced by the solid rheology, while with relatively low fluid mobility (K00 ∼10−11 m3 s kg−1) and thus a low compaction length (δ0 ∼ 2 km), the influence of solid advection becomes important.

In our recent work (Cerpa et al., 2017), we have investigated the combined effects of grain-size distribution and mantle wedge viscosity structure on fluid-flow pathways. In these models, the reference mobility is ∼10–9 m3 s kg−1, which is equivalent to moderate mobility in Wilson et al. (2014). However, the spatial variation in grain size caused the fluid mobility and thus the compaction length to vary spatially. Fluid mobility is relatively low near the MDD where grain size is small and increases with the distance away from the MDD. As a result, fluids introduced beneath the forearc are initially trapped at the base of the mantle wedge and dragged downdip by mantle flow. The downdip increase in permeability increases fluid mobility and thus decreases the effectiveness of solid advection, allowing the upward migration of fluids at subarc depths as shown in our model F1 (Fig. 4). In other regions of the mantle wedge, where permeability is relatively high, fluids are largely unaffected by solid advection except in regions where variations in solid shear viscosity promote lateral fluid migration through its effects on compaction pressure gradients and solid advection. In the latter regions, some fluids migrate toward the trench and may contribute to the focusing of some fluids released at postarc depths toward the subarc region. Overall, spatial variations in both grain size and mantle viscosity help to focus fluids beneath the arc, but their contribution to the focusing varies spatially.

In the present work, we have investigated the effects of fluid influx in terms of ω and fluid viscosity on fluid migration pathways. The results illustrate the scaling relationship between ω and µ such that the fluid migration pathways are identical in models with the same ω2/µ value, although the distribution of fluid fraction varies with ω. A low ω2/µ value (2.25 × 10−6 Pa−1 s−1) leads to the entrapment of fluids at the base of the mantle wedge and their downdip transport by the mantle flow, inconsistent with the predictions of flux melting occurring beneath the arcs. A very high ω2/µ value (2.5) leads to mostly upward fluid migration. At intermediate to high ω2/µ values (2.5 × 10−4–2.5 × 10−2 Pa−1 s−1), fluids released beneath the forearc are transported downdip to the subarc depths by solid advection before migrating upward, and some fluids released during a second peak dehydration at postarc depths are laterally advected toward the arc in the mantle inflow region. Increasing ω2/µ (>2.5 × 10−2 Pa−1 s−1) reduces the lateral transport of fluids by solid advection and the widening of the pathways of rising fluids through the mantle wedge beneath the backarc.

Wilson et al. (2014) observed the formation of updip fluid-flow channels in the subducting slab in models with high fluid mobility and thus a high compaction length. Our models predict, instead, that some fluids released at postarc depths can flow updip near the base of the mantle wedge if ω2/µ is relatively high, providing a new mechanism for focusing fluids beneath the subarc region. Geophysical observations in NE Japan suggest the existence of a dipping low-velocity layer near the base of the mantle wedge (Kawakatsu and Watada, 2007. This layer may indicate the presence of free fluids, but whether fluids are migrating downdip or updip would be difficult to resolve. Our models show that the direction of fluid migration at the base of the mantle wedge at postarc depths depends on ω2/µ. Alternatively, it has been proposed that the low-velocity layer just above the slab represents the downdip drag of hydrous phases that are stable up to 150 km depth (Iwamori, 1998; Grove et al., 2009; Horiuchi and Iwamori, 2016), although thermomechanical models with realistic rheologies for the mantle wedge predict that the temperatures are too high for most hydrous phases to remain stable beyond a depth of 70–80 km (e.g., van Keken et al., 2011; Wada et al., 2012).

A common simplification made in fluid migration models for subduction zone is the use of a constant fluid viscosity as in our models F1 to F7. In nature, however, fluid viscosity may change spatially because of its dependence on pressure, temperature, and composition of the fluids. As a first step in incorporating fluid viscosity variation, we applied temperature-dependent fluid viscosity in models F8 and F9. The results show that the relatively cold temperature just above the slab and the subsequent low mobility may contribute to the drag of fluid released beneath the forearc toward the subarc region. However, the fluid migration paths in these models do not differ significantly from those in the models with a constant fluid viscosity. Production of silicic melts in the fluxed region may produce important changes in fluid properties but is neglected in our models. Hence, fluid migration pathways in the shallow mantle wedge may differ from the predictions by our models.

Although a relatively large effect was observed when the fluid density was decreased by one order of magnitude, it is unlikely to decrease by more than a factor of two in reality (Hack and Thompson, 2011). We can thus conclude that in a realistic range of fluid density in the mantle wedge, its effect on fluid migration is small compared to the effects of fluid viscosity and fluid influx on fluid migration pathways. Nonetheless, the results suggest that the combination of high fluid density and fluid viscosity due to high silicate content in the fluxed region of the mantle can have a large impact on fluid flow in the shallow part of the mantle wedge, enhancing the effect of solid advection. Hence, as we have suggested in Cerpa et al. (2017), the focusing of fluids toward the subarc region may be greater than what is predicted by our models.

6.2. Fluid Distribution and Hydrous Melting

The supply of water from the subducting slab to the hot mantle wedge is thought to trigger hydrous melting (Sisson and Grove, 1993; Grove et al., 2006), and the degree of melting likely increases with the amount of water added to the system (Stolper and Newman, 1994; Hirschmann et al., 1999; Katz et al., 2003; Kelley et al., 2010). Although, our models do not incorporate melting in the calculations, we use the parameterized melting function of Kelley et al. (2010), following the approach of Wilson et al. (2014), to calculate the fraction of hydrous melts that can be produced in the mantle wedge based on the results of selected fluid-flow models that generate some upward fluid migration and time-averaged porosities ≥10−4 (models F1, F2, F3, F6, and F7). The degree of melting is calculated by using fluid content, temperature, and pressure predicted by our fluid migration models. The parameterized melting function is given in terms of temperature, pressure, and water content. We use the temperature field from model S1 (Fig. 2). A lithostatic pressure based on solid density of 3300 kg m−3 (Δρ = 2000 kg m−3) is assumed for the calculation of pressure. Water concentration is calculated using the time-averaged fluid fraction of the fluid-flow models. We exclude from the calculations the shallowest part of the mantle wedge where fluid accumulation due to ponding is likely overestimated by our models without inclusion of the effects of brittle deformation on fluid extraction (Cerpa et al., 2017).

The fluid-flow models predict the formation of two or three main pathways through the hottest part of the mantle wedge. The first main fluid pathway forms at relatively shallow depth (<100 km) beneath the arc (models F1 and F6) or beneath the forearc at a distance up to ∼25 km from the arc, if fluid mobility is high (model F3). A second main fluid pathway is observed in all models beneath the backarc region and generally forms above the location of the second peak dehydration situated between 135 and 150 km depths. A third main pathway eventually appears beneath the backarc region as a consequence of the updip flow of some fluids from the second peak dehydration, if fluid mobility is high (models F2, F3, and F7).

Our calculations predict that hydrous melting occurs beneath the arc and the backarc (Fig. 12). Beneath the arc, melt production is observed in a relatively narrow region (width <20 km) between 70 and 85 km depths near the updip sloped isocontours of solid shear viscosity in the inflow region of the mantle wedge. Beneath the backarc, melt production develops in a wider region than in the subarc as the fluid-flow calculations predict wider fluid pathways in the former region. Overall, melt fraction is relatively low (0.1%–1%), reaching 1%–4% only in models F5 and F7. At a given value of ω2/µ (identical fluid pathways), the highest melt fractions are observed for the models with the highest ω.

The degree of melting in the mantle wedge inferred from geochemical data on arc lavas and volcanics ranges from 1% to 20% (Portnyagin et al., 2007; Kelley et al., 2010; Watt et al., 2013). Our models generally predict a lower degree of melting in the mantle wedge (≤1%). This can be explained by underestimation of the amount of slab-derived fluids in the hottest part of the mantle wedge or the temperature or both. Increasing tuning parameter increases the predicted melt fraction but requires relatively high fluid viscosity to cause focusing of fluids beneath the arc (model F7). Alternatively, mantle potential temperature may be higher than that assumed in models F1–F11 (Tm = 1350 °C) and may be closer to the dry solidus temperature of mantle peridotite. Moreover, experimental studies on melting of natural peridotite samples in the presence of water indicate that low-degree melts may be produced at the base of the mantle wedge at relatively low temperatures of ∼800–900 °C below an 80 km depth (Grove et al., 2006; Till et al., 2012). The latent heat released by such melting would be advected by the rising melts, locally increasing the temperature in the core part of the mantle wedge and promoting further melting (Rees-Jones et al., 2018).

Our models with relatively old subducting slabs also predict hydrous melting beneath the backarc. In nature, volcanoes do form behind the volcanic front above relatively old subducting slabs, such as in Northeast Japan (Tatsumi et al., 1983), Marianas (Kelley et al. 2010), and Kamchatka (e.g., Tatsumi et al., 1994), indicating the presence of magma as predicted by our models As discussed in Cerpa et al. (2017), some of the melts that are produced beneath the backarc may also be redirected arcward by the incoming mantle flow, particularly when the fluid viscosity and density are relatively high (i.e., low fluid mobility and buoyancy). Gradients of compaction pressure induced by melting alone may also facilitate lateral transport of melts as proposed for melt focusing at mid-ocean ridges (Turner et al., 2017).

6.3. Northern Cascadia

Our fluid-flow models for Cascadia (models FC11–FC13) (Fig. 11) display fluid pathways comparable to those inferred from MT models. Our solid-state flow model for Cascadia predicts the release of fluids into the flowing mantle wedge between 75 and 95 km depths, consistent with the highest electrical conductivity observed just above the slab between 80 and 100 km depths by (McGary et al., 2014). The high-conductivity anomaly at the base of the mantle wedge extends downdip up to a 120 km depth, and there appears to be little fluid present beyond the 120 km depth, if any. Similarly, our models with the highest fluid mobility (FC11 and FC13) exhibit modest fluid fraction beyond 110 km depth. However, the MT model of Wannamaker et al. (2014) for the same Cascadia profile indicates another slab-parallel high-conductivity anomaly at the base of the mantle wedge between a 150 km and a 200 km depth. These authors interpreted this anomaly as a region of decompression melting of upwelling mantle. Our models (models FC11–FC13) predict that some fluids do become dragged downdip by mantle flow to greater depth, potentially contributing to the deeper high-conductivity anomaly.

In all three Cascadia models, despite two different fluid mobilities (FC11 and FC13 versus FC12), upward migration of fluids occurs in the hottest part of the mantle wedge near the subarc region as observed by McGary et al. (2014) and Wannamaker et al. (2014). Their MT models indicate some lateral deflection of upward fluid pathways that may be explained by the combined effect of mantle inflow and vertical variations in mantle shear viscosity as shown in our models. However, our Cascadia models display upward fluid pathways east of Mount Rainier, while the MT models indicate upward pathways west of Mount Rainier. This discrepancy points to two possible explanations. One is that the depth of fluid influx is shallower in the actual system, potentially due to higher-than-expected slab temperature and thus shallower fluid release or the updip fluid migration within the subducting slab that leads to fluid influx at shallower depths than in the model. The other possible explanation is that the permeability structure at the base of the mantle wedge updip of the influx region is higher in the actual system, for example, due to shear-induced permeability (e.g., Précigout et al., 2017) or reaction-induced fracturing (Okamoto and Shimizu, 2015), allowing updip migration along the base of the mantle wedge.

We use the time-averaged porosity distribution of our Cascadia models to predict the regions of hydrous melting (Fig. 13). The model with the lowest fluid viscosity (model FC11) predicts negligible melting. The two models with higher fluid viscosity (models FC12 and FC13), predict 0.1%–1% degree of melting in the subarc region. This relatively low degree of melting predicted by these models is consistent with the relatively high ratios of incompatible to compatible elements found in Cascadia arc lavas (e.g., Ce/Yb; Schmidt et al., 2008). In hot-subduction zones such as Cascadia, relatively high fluid viscosity is expected because of the temperature dependence of silicates solubility in fluids. Further, fluids released from the dehydration of the subducting mantle in the subarc region may trigger slab melting as they travel through the overlying crust as proposed to explain geochemical data in Southern Cascadia (Walowski et al., 2015), increasing the viscosity of the fluids before they migrate into the mantle wedge. Thus, as indicated by models FC12 and FC13 (Fig. 13), with relatively high fluid viscosity, a modest fraction of slab-derived fluids from the slab can explain 0.1%–1% melting beneath Mount Rainier, and melting occurs in a relatively small region compared to the distribution of fluids in the mantle wedge (Fig. 11).

7. CONCLUSIONS

In our previous study (Cerpa et al. [2017]), we showed that the spatial variations in grain size (permeability) and mantle shear viscosity play a key role on fluid migration in the mantle wedge. In the present study, we assess the influence of the magnitude of fluid influx through the parameter ω and fluid properties in our fluid migration models. In these models, fluids are generally introduced from the base of the mantle wedge beneath the forearc and the backarc. Our results suggest that fluid influx and fluid viscosity (µ) have a large effect on fluid migration pathways, while density contrast plays a secondary role. In our models, mantle compaction is mainly controlled by ω2/µ for a given grain-size distribution. Low values of ω2/µ produce entrapment of fluids at the base of the mantle wedge. Increasing ω2/µ allows upward migration of fluids. In particular, the downdip drag of fluid released beneath the forearc decreases with increasing ω2/µ. At high values, the models predict updip flow of fluids at the base of the mantle wedge beneath the backarc. Further, temperature-dependent fluid viscosity enhances the downdip drag of fluids by mantle flow. Based on the averaged porosities predicted by our models, we estimate relatively low degree of melting (0.1%–5%) occurring beneath the arc and beneath the backarc. Finally, our model for Northern Cascadia predicts fluid distribution in the mantle wedge that is consistent with electromagnetic imaging of the mantle wedge.

ACKNOWLEDGMENTS

We thank T. Gerya and an anonymous reviewer for comments that helped to improve the manuscript. I.W. acknowledges financial support from the University of Minnesota in the form of start-up funds, which also supported N.G.C. C.W. was supported by NSF grant OCE-1358091. The code TerraFERMA is open-source and available at http://terraferma.github.io.

APPENDIX A. PARAMETERS IN SOLID-FLOW CALCULATIONS

Solid Shear Viscosity

The solid-state mantle flow model accounts for diffusion creep and dislocation creep, which are assumed to be independent. The solid shear viscosity η is thus expressed as 
graphic
where ηmax is a maximum shear viscosity set to 1024 Pa s, and ηdiff and ηdisl are diffusion and dislocation shear viscosity, respectively, and are defined as 
graphic
where the subscript i denotes creep mechanism, b is grain size, and ⋵II is the second invariant of strain rate. The rheological parameters Λ, COH, E, and R are pre-exponential factor, water content, activation energy, and gas constant, respectively. Parameters α, β, and γ are grain size exponent, water content exponent, and stress exponent, respectively. For these parameters, we use values that are determined experimentally for wet olivine (Hirth and Kohlstedt, 2003) (Table A1).

Grain Size

We calculate the steady-state grain size as (Behn et al., 2009) 
graphic
where σ and ⋵disl are stress tensor and strain rate tensor associated with dislocation creep, respectively, G is grain growth constant, λ is average specific grain boundary energy, Eg is activation energy for grain growth, C is a geometric constant, pg is grain growth exponent, and χ represents the fraction of work done by dislocation creep associated with changing grain boundary area. The value of the parameters used for the calculations are provided Table A2.

APPENDIX B. SCALING RELATIONSHIPS

The reference compaction length is given by Spiegelman (1993): 
graphic
where we have defined the reference porosity as: 
graphic
with Ψ0 a reference flux used for the non-dimensionalization of the boundary condition at the base of the mantle wedge following: 
graphic
Note that the reference flux and density contrast are also used to define the reference fluid velocity: 
graphic

The reference permeability (K0 = 3.7 × 10−9m2), fluid viscosity (µ0 = 1 Pa s), and solid viscosity (η0 = 1019 Pa s) are fixed in all models.

1Supplemental Material. Data used to generate figures. Please visit https://doi.org/10.1130/GES01660.S1 or access the full-text article on www.gsapubs.org to view the supplemental data. The code TerraFERMA is open-source and available at http://terraferma.github.io.
Science Editor: Shanaka de Silva
Associate Editor: Philippe Agard
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.